Added support for static libraries,
[u/mrichter/AliRoot.git] / HLT / trigger / RunD0Trigger.C
1 //#define __COMPILE__
2 //#ifdef __COMPILE__
3 #if !defined(__CINT__) || defined(__MAKECINT__)
4 //-- --- standard headers------------- 
5 #include <iostream.h>
6 //--------Root headers ---------------
7 #include <TSystem.h>
8 #include <TFile.h>
9 #include <TString.h>
10 #include <TStopwatch.h>
11 #include <TObject.h>
12 #include <TVector3.h>
13 #include <TTree.h>
14 #include <TParticle.h>
15 #include <TArray.h>
16 //----- AliRoot headers ---------------
17 #include "alles.h"
18 #include "AliRun.h"
19 #include "AliKalmanTrack.h"
20 #include "AliITStrackV2.h"
21 #include "AliHeader.h"
22 #include "AliGenEventHeader.h"
23 #include "AliV0vertex.h"
24 #include "AliV0vertexer.h"
25 #include "AliITSVertex.h"
26 #include "AliITSVertexer.h"
27 #include "AliITSVertexerTracks.h"
28 #include "AliD0Trigger.h"
29 #endif
30 //-------------------------------------
31 // field (T)
32 const Double_t kBz = 0.4;
33
34 // primary vertex
35 Double_t primaryvertex[3] = {0.,0.,0,};
36
37 // sigle track cuts
38 const Double_t kPtCut = 0.5;  // GeV/c
39 const Double_t kd0Cut = 50.; // micron
40 const Double_t kd0CutHigh = 200.; // micron
41
42
43 //cuts for combined tracks
44 // cuts[0] = lowest V0 cut
45 // cuts[1] = highest V0 cut
46 // cuts[2] = inv. mass cut (diiferense)
47 // cuts[3] = max value for pointing angle
48 const Double_t cuts[4] = {0.005,
49                           0.02,
50                           0.5,
51                           0.95};
52
53 // this function applies single track cuts
54 Bool_t TrkCuts(const AliITStrackV2& trk);
55
56 // this function creates TObjArrays with positive and negative tracks
57 void   SelectTracks(TTree& itsTree,
58                       TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
59                       TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
60
61 //void GetPrimaryVertex(int i,Char_t* path="./");
62
63 Int_t iTrkP,iTrkN,itsEntries;
64 Char_t trksName[100];
65 Int_t nTrksP=0,nTrksN=0;
66 Int_t nD0=0;
67 int ev=0;
68
69 void   RunD0Trigger(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
70
71   const Char_t *name="AliD0Trigger";
72   cerr<<'\n'<<name<<"...\n";
73   gBenchmark->Start(name); 
74
75   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
76
77   // Open file with ITS tracks
78   Char_t fnameTrack[1024];
79   sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
80   TFile* itstrks = TFile::Open(fnameTrack);
81
82    // tracks from ITS
83   sprintf(trksName,"TreeT_ITS_%d",ev);
84   TTree *itsTree=(TTree*)itstrks->Get(trksName);
85   if(!itsTree) continue;
86   itsEntries = (Int_t)itsTree->GetEntries();
87   printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
88
89   //Getting primary Vertex
90   GetPrimaryVertex(0,path);
91
92   // call function which applies sigle track selection and
93   // separetes positives and negatives
94   TObjArray trksP(itsEntries/2);
95   Int_t *itsEntryP = new Int_t[itsEntries];
96   TObjArray trksN(itsEntries/2);
97   Int_t *itsEntryN = new Int_t[itsEntries];
98   SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN); 
99
100   cout<<"#pos: "<<nTrksP<<endl;
101   cout<<"#neg: "<<nTrksN<<endl;
102
103   AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
104
105   for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
106     postrack = (AliITStrackV2*)trksP.At(iTrkP);
107     for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
108       negtrack = (AliITStrackV2*)trksN.At(iTrkN);
109       D0.SetTracks(postrack,negtrack);
110       if(D0.FindV0()){
111         D0.FindMomentaAtVertex();
112         if(D0.FindInvMass()){
113           if(D0.PointingAngle()){
114             nD0++;
115           }
116         }
117       } 
118     }
119   }
120   cout<<"#D0: "<<nD0<<endl;
121   gBenchmark->Stop(name);
122   gBenchmark->Show(name);
123 }
124 //___________________________________________________________________________
125 void   SelectTracks(TTree& itsTree,
126                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
127                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
128 //
129 // this function creates two TObjArrays with positive and negative tracks
130 //
131   nTrksP=0,nTrksN=0;
132
133  
134   Int_t entr = (Int_t)itsTree.GetEntries();
135
136   // trasfer tracks from tree to arrays
137   for(Int_t i=0; i<entr; i++) {
138
139     AliITStrackV2 *itstrack = new AliITStrackV2; 
140     itsTree.SetBranchAddress("tracks",&itstrack);
141
142     itsTree.GetEvent(i);
143
144     // single track selection
145     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
146
147     if(itstrack->Get1Pt()>0.) { // negative track
148       trksN.AddLast(itstrack);
149       itsEntryN[nTrksN] = i;
150       nTrksN++;
151     } else {                    // positive track
152       trksP.AddLast(itstrack);
153       itsEntryP[nTrksP] = i;
154       nTrksP++;
155     }
156
157   }
158
159   return;
160 }
161 //____________________________________________________________________________
162 Bool_t TrkCuts(const AliITStrackV2& trk) {
163 // 
164 // this function tells if track passes some kinematical cuts  
165 //
166   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
167   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
168   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
169
170   return kTRUE;
171 }
172 //____________________________________________________________________________
173 void GetPrimaryVertex(int i,Char_t* path="./") {
174
175   int event=i;
176
177   Char_t falice[1024];
178   sprintf(falice,"%s/galice.root",path);
179   TFile * galice = new TFile(falice);
180   
181   TDirectory * curdir;  
182
183   Char_t vname[20];
184   galice->cd();
185   
186   sprintf(vname,"Vertex_%d",event);
187   TArrayF o = 0;
188   o.Set(3);
189   AliHeader * header = 0;
190   
191   TTree * treeE = (TTree*)gDirectory->Get("TE");
192   treeE->SetBranchAddress("Header",&header);
193   treeE->GetEntry(event);
194   AliGenEventHeader* genHeader = header->GenEventHeader();
195   if(genHeader){
196     genHeader->PrimaryVertex(o);
197     primaryvertex[0] = (Double_t)o[0];
198     primaryvertex[1] = (Double_t)o[1];
199     primaryvertex[2] = (Double_t)o[2];
200   }
201   else{
202     printf("Can't find Header");
203   }
204   delete header;
205   delete galice;
206 }