Merged Bergen, mergen Cvetan TransformerRow and
[u/mrichter/AliRoot.git] / HLT / trigger / RunD0offline.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 primaryvertex[3]={0.,0.,0.};
36
37 //sec. vertex
38 double v2[3]={0,0,0};
39
40 // sigle track cuts
41 const Double_t kPtCut = 0.5;        // 0.5 GeV/c
42 const Double_t kd0Cut = 50.;        // 50  micron
43 const Double_t kd0CutHigh = 400.;   // 200 micron
44
45
46 //cuts for combined tracks
47 const Double_t cuts[7] = {0.005,     // 0.005 cuts[0] = lowest V0 cut  (cm)
48                           0.800,     // 0.015 cuts[1] = highest V0 cut (cm)
49                           0.012,     // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
50                           0.8,       // 0.8   cuts[3] = min. cosine value for pointing angle
51                           -5000,     // -5000 cuts[4] = d0d0
52                           0,       // 0.8   cuts[5] = costhetastar
53                           0.5};      // 0.5   cuts[6] = ptchild
54 //cut for distance of closest aprach
55 double cutDCA=0.05;   //0.05
56
57 // this function applies single track cuts
58 Bool_t TrkCuts(const AliITStrackV2& trk);
59
60 // this function creates TObjArrays with positive and negative tracks
61 void   SelectTracks(TTree& itsTree,
62                       TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
63                       TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
64
65 //void GetPrimaryVertex(int i,Char_t* path="./");
66
67 //void PtD0(Char_t* path="./");
68
69 Int_t iTrkP,iTrkN,itsEntries;
70 Char_t trksName[100];
71 Int_t nTrksP=0,nTrksN=0;
72 Int_t nD0=0;
73 int ev=0;
74 double mom[6];
75 int event[10000];
76 int index=0;
77
78 void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
79   Char_t falice[1024];
80   sprintf(falice,"%s/galice.root",path);
81   TFile *f = new TFile(falice);   
82   gAlice=(AliRun*)f->Get("gAlice");
83   int nEvent=gAlice->GetEventsPerRun();
84   //int nEvent=20;
85   delete gAlice;
86
87   const Char_t *name="AliD0offline";
88   cerr<<'\n'<<name<<"...\n";
89   cout<<"#Events: "<<nEvent<<endl;
90   gBenchmark->Start(name); 
91
92   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
93
94   // Open file with ITS tracks
95   Char_t fnameTrack[1024];
96   //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
97   sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
98   TFile* itstrks = TFile::Open(fnameTrack);
99
100   //the offline stuff
101   // define the cuts for vertexing
102   Double_t vtxcuts[]={33., // max. allowed chi2
103                       0.0, // min. allowed negative daughter's impact param 
104                       0.0, // min. allowed positive daughter's impact param 
105                       1.0, // max. allowed DCA between the daughter tracks
106                       -1.0, // min. allowed cosine of V0's pointing angle
107                       0.0, // min. radius of the fiducial volume
108                       2.9};// max. radius of the fiducial volume
109
110   TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
111   TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with |Eta|<0.9",100,0,10);
112   TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
113   TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
114
115   for(ev=0;ev<nEvent;ev++) {
116     
117     cout<<"\n Event: "<<ev<<endl;
118
119     // tracks from ITS
120     sprintf(trksName,"TreeT_ITS_%d",ev);
121     TTree *itsTree=(TTree*)itstrks->Get(trksName);
122     if(!itsTree) continue;
123     itsEntries = (Int_t)itsTree->GetEntries();
124     printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
125     
126     //Getting primary Vertex
127     GetPrimaryVertex(ev,path);
128     
129     // call function which applies sigle track selection and
130     // separetes positives and negatives
131     TObjArray trksP(itsEntries/2);
132     Int_t *itsEntryP = new Int_t[itsEntries];
133     TObjArray trksN(itsEntries/2);
134     Int_t *itsEntryN = new Int_t[itsEntries];
135     SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN); 
136     
137     cout<<"#pos: "<<nTrksP<<endl;
138     cout<<"#neg: "<<nTrksN<<endl;
139     
140     // create the AliV0vertexer object
141     AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
142     
143     AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
144     
145     double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
146     
147     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
148       postrack = (AliITStrackV2*)trksP.At(iTrkP);
149       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
150         negtrack = (AliITStrackV2*)trksN.At(iTrkN);
151         D0->SetTracks(postrack,negtrack);
152         //
153         // ----------- DCA MINIMIZATION ------------------
154         //
155         // find the DCA and propagate the tracks to the DCA 
156         double dca = vertexer2->PropagateToDCA(negtrack,postrack);
157         
158         if(dca<cutDCA){
159           // define the AliV0vertex object
160           AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
161           // get position of the vertex
162           vertex2->GetXYZ(v2[0],v2[1],v2[2]);
163           delete vertex2;       
164           //if(D0->FindV0offline(v2) && D0->d0d0()){
165           if(D0->d0d0()){
166             D0->SetV0(v2);
167             D0->FindMomentaOffline();
168             if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
169               nD0++;
170               event[index]=ev; index++;
171               if(h){
172                 h1->Fill(D0->Pt());
173                 h3->Fill(D0->Eta());
174               }
175             }
176           } 
177         }
178       }
179     } 
180   }
181   cout<<"\n#D0: "<<nD0<<endl;
182   for(int i=0;i<nD0;i++)
183     cout<<event[i]<<endl;
184   gBenchmark->Stop(name);
185   gBenchmark->Show(name);
186
187   if(h){
188     if(!PtD0){
189       TCanvas *c = new TCanvas("c","c",700,1000);
190       c->Divide(1,2);
191       c->cd(1);
192       h1->Draw();
193       pt = new TPaveText(7.3,0.4,11,1.8,"br");
194       pt->SetFillColor(0);
195       pt->SetBorderSize(1);
196       pt->AddText("Cuts:");
197       Char_t st[1024];
198       sprintf(st,"First Pt:         %g",kPtCut);
199       pt->AddText(st);
200       sprintf(st,"d0 low:           %g",kd0Cut);
201       pt->AddText(st);
202       sprintf(st,"d0 high:          %g",kd0CutHigh);
203       pt->AddText(st);
204       sprintf(st,"V0 low:           %g",cuts[0]);
205       pt->AddText(st);
206       sprintf(st,"V0 high:          %g",cuts[1]);
207       pt->AddText(st);
208       sprintf(st,"InvMass Diff:     %g",cuts[2]);
209       pt->AddText(st);
210       sprintf(st,"cosPointingAngle: %g",cuts[3]);
211       pt->AddText(st);
212       sprintf(st,"d0d0:             %g",cuts[4]);
213       pt->AddText(st);
214       sprintf(st,"cosThetaStar:     %g",cuts[5]);
215       pt->AddText(st);
216       sprintf(st,"PtChild:          %g",cuts[6]);
217       pt->AddText(st);
218       sprintf(st,"DCA:              %g",cutDCA);
219       pt->AddText(st);
220       pt->Draw();
221       c_1->Modified();
222       c->cd();
223       c->cd(2);
224       h3->Draw();
225     }
226     if(PtD0){
227       PtD0(path,h2,h4);
228       TCanvas *c = new TCanvas("c","c",1000,700);
229       c->Divide(2,2);
230       c->cd(1);
231       h1->Draw();
232       pt = new TPaveText(7.3,0.4,11,1.8,"br");
233       pt->SetFillColor(0);
234       pt->SetBorderSize(1);
235       pt->AddText("Cuts:");
236       Char_t st[1024];
237       sprintf(st,"First Pt:         %g",kPtCut);
238       pt->AddText(st);
239       sprintf(st,"d0 low:           %g",kd0Cut);
240       pt->AddText(st);
241       sprintf(st,"d0 high:          %g",kd0CutHigh);
242       pt->AddText(st);
243       sprintf(st,"V0 low:           %g",cuts[0]);
244       pt->AddText(st);
245       sprintf(st,"V0 high:          %g",cuts[1]);
246       pt->AddText(st);
247       sprintf(st,"InvMass Diff:     %g",cuts[2]);
248       pt->AddText(st);
249       sprintf(st,"cosPointingAngle: %g",cuts[3]);
250       pt->AddText(st);
251       sprintf(st,"d0d0:             %g",cuts[4]);
252       pt->AddText(st);
253       sprintf(st,"cosThetaStar:     %g",cuts[5]);
254       pt->AddText(st);
255       sprintf(st,"PtChild:          %g",cuts[6]);
256       pt->AddText(st);
257       sprintf(st,"DCA:              %g",cutDCA);
258       pt->AddText(st);
259       pt->Draw();
260       c_1->Modified();
261       c->cd();
262       c->cd(2);
263       h3->Draw();
264       c->cd(3);
265       h2->Draw();
266       c->cd(4);
267       h4->Draw();
268     }
269     //TFile *outfile = new TFile("results.root","recreate");
270     //h1->Write();
271     //outfile->Close();
272     
273   } 
274 }
275 //___________________________________________________________________________
276 void   SelectTracks(TTree& itsTree,
277                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
278                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
279 //
280 // this function creates two TObjArrays with positive and negative tracks
281 //
282   nTrksP=0,nTrksN=0;
283
284  
285   Int_t entr = (Int_t)itsTree.GetEntries();
286
287   // trasfer tracks from tree to arrays
288   for(Int_t i=0; i<entr; i++) {
289
290     AliITStrackV2 *itstrack = new AliITStrackV2; 
291     itsTree.SetBranchAddress("tracks",&itstrack);
292
293     itsTree.GetEvent(i);
294
295     // single track selection
296     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
297
298     if(itstrack->Get1Pt()>0.) { // negative track
299       trksN.AddLast(itstrack);
300       itsEntryN[nTrksN] = i;
301       nTrksN++;
302     } else {                    // positive track
303       trksP.AddLast(itstrack);
304       itsEntryP[nTrksP] = i;
305       nTrksP++;
306     }
307
308   }
309
310   return;
311 }
312 //____________________________________________________________________________
313 Bool_t TrkCuts(const AliITStrackV2& trk) {
314 // 
315 // this function tells if track passes some kinematical cuts  
316 //
317   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
318   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
319   //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
320
321   return kTRUE;
322 }
323 //____________________________________________________________________________
324 void GetPrimaryVertex(int i,Char_t* path="./") {
325
326   int event=i;
327
328   Char_t falice[1024];
329   sprintf(falice,"%s/galice.root",path);
330   TFile * galice = new TFile(falice);
331   
332   TDirectory * curdir;  
333
334   Char_t vname[20];
335   galice->cd();
336   
337   sprintf(vname,"Vertex_%d",event);
338   TArrayF o = 0;
339   o.Set(3);
340   AliHeader * header = 0;
341   
342   TTree * treeE = (TTree*)gDirectory->Get("TE");
343   treeE->SetBranchAddress("Header",&header);
344   treeE->GetEntry(event);
345   AliGenEventHeader* genHeader = header->GenEventHeader();
346   if(genHeader){
347     genHeader->PrimaryVertex(o);
348     primaryvertex[0] = (Double_t)o[0];
349     primaryvertex[1] = (Double_t)o[1];
350     primaryvertex[2] = (Double_t)o[2];
351   }
352   else{
353     printf("Can't find Header");
354   }
355   delete header;
356   delete galice;
357 }
358 //____________________________________________________________________________
359 PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
360   
361   Char_t falice[1024];
362   sprintf(falice,"%s/galice.root",path);  
363   TFile *f = new TFile(falice);
364   gAlice=(AliRun*)f->Get("gAlice");
365   
366   TParticle *p;
367   int nd0=0;
368   int nkminus =0;
369   int npipluss = 0;
370   int nEvent=gAlice->GetEventsPerRun();
371    
372   for (Int_t i = 0; i <nEvent; i++) {
373     cout<<"Event: "<<i<<endl;
374     gAlice->GetEvent(i);
375     Int_t nPart = gAlice->GetNtrack();
376     for (Int_t iPart = 0; iPart < nPart; iPart++) {
377       //cout<<"Particlenr.: "<<iPart<<endl;
378       p = (TParticle*)gAlice->Particle(iPart);
379       if (p->GetPdgCode()==421){
380         if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
381         h4->Fill(p.Eta());
382         nd0++;
383       }
384       if (p->GetPdgCode()==-321){
385         nkminus++;
386       }
387       if (p->GetPdgCode()==211){
388         npipluss++;
389       }
390     }
391   }
392   delete gAlice;
393 }