Updated macros for PHOS alignment calculation
[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 #include <fstream.h>
7 //--------Root headers ---------------
8 #include <TSystem.h>
9 #include <TFile.h>
10 #include <TString.h>
11 #include <TStopwatch.h>
12 #include <TObject.h>
13 #include <TVector3.h>
14 #include <TTree.h>
15 #include <TObjArray>
16 #include <TParticle.h>
17 #include <TArray.h>
18 //----- AliRoot headers ---------------
19 #include "alles.h"
20 #include "AliRun.h"
21 #include "AliKalmanTrack.h"
22 #include "AliITStrackV2.h"
23 #include "AliHeader.h"
24 #include "AliGenEventHeader.h"
25 #include "AliV0vertex.h"
26 #include "AliV0vertexer.h"
27 #include "AliITSVertex.h"
28 #include "AliITSVertexer.h"
29 #include "AliITSVertexerTracks.h"
30 #include "AliD0Trigger.h"
31 #endif
32 //-------------------------------------
33 typedef struct {
34   Int_t lab;
35   Int_t pdg;
36   Int_t mumlab;
37   Int_t mumpdg;
38   Float_t Vx,Vy,Vz;
39   Float_t Px,Py,Pz;
40 } RECTRACK;
41
42 // field (T)
43 const Double_t kBz = 0.4;
44
45 // primary vertex
46 double primaryvertex[3]={0.,0.,0.};
47
48 //sec. vertex
49 double v2[3]={0,0,0};
50
51 // sigle track cuts
52 const Double_t kPtCut = 0.5;        // 0.5 GeV/c
53 const Double_t kd0Cut = 50.;        // 50  micron
54 const Double_t kd0CutHigh = 400.;   // 200 micron
55
56
57 //cuts for combined tracks
58 const Double_t cuts[7] = {0.005,     // 0.005 cuts[0] = lowest V0 cut  (cm)
59                           0.800,     // 0.015 cuts[1] = highest V0 cut (cm)
60                           0.012,     // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
61                           0.8,       // 0.8   cuts[3] = min. cosine value for pointing angle
62                           -5000,     // -5000 cuts[4] = d0d0
63                           0,         // 0.8   cuts[5] = costhetastar
64                           0.5};      // 0.5   cuts[6] = ptchild
65 //cut for distance of closest aprach
66 double cutDCA=0.05;   //0.05
67
68 // this function applies single track cuts
69 Bool_t TrkCuts(const AliITStrackV2& trk);
70
71 // this function creates TObjArrays with positive and negative tracks
72 void   SelectTracks(TTree& itsTree,
73                       TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
74                       TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
75
76 //void GetPrimaryVertex(int i,Char_t* path="./");
77
78 //void PtD0(Char_t* path="./");
79
80 Int_t iTrkP,iTrkN,itsEntries;
81 Char_t trksName[100],refsName[100];
82 Int_t nTrksP=0,nTrksN=0;
83 Int_t nD0=0;
84 int ev=0;
85 double mom[6];
86 int event[10000];
87 int index=0;
88 Bool_t isSignal;
89 Int_t nTotEv=0,nD0recSgn=0,nD0recBkgS=0,nD0recBkg=0,nD0rec1ev=0;
90 Int_t pdg[2],mum[2],mumlab[2];
91 RECTRACK rectrk;
92
93 void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
94     
95   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
96   
97   // Open file with ITS tracks
98   Char_t fnameTrack[1024];
99   //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
100   sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
101   TFile* itstrks = TFile::Open(fnameTrack);
102   
103   Char_t refFile[1024];
104   //sprintf(fnameTrack,"%s/ITStracksRefFile.root",path);
105   sprintf(refFile,"%s/ITStracksRefFile.root",path);
106   TFile* itsrefs = TFile::Open(refFile);
107   
108
109   //the offline stuff
110   // define the cuts for vertexing
111   Double_t vtxcuts[]={33., // max. allowed chi2
112                       0.0, // min. allowed negative daughter's impact param 
113                       0.0, // min. allowed positive daughter's impact param 
114                       1.0, // max. allowed DCA between the daughter tracks
115                       -1.0, // min. allowed cosine of V0's pointing angle
116                       0.0, // min. radius of the fiducial volume
117                       2.9};// max. radius of the fiducial volume
118   
119   TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
120   TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with |Eta|<0.9",100,0,10);
121   TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
122   TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
123   
124   Char_t falice[1024];
125   sprintf(falice,"%s/galice.root",path);
126   TFile *f = new TFile(falice);   
127   gAlice=(AliRun*)f->Get("gAlice");
128   int nEvent=gAlice->GetEventsPerRun();
129   //int nEvent=20;
130   cout<<"#Events: "<<nEvent<<endl;
131   delete gAlice;
132
133   const Char_t *name="AliD0offline";
134   cerr<<'\n'<<name<<"...\n";
135   gBenchmark->Start(name); 
136   
137   for(ev=0;ev<nEvent;ev++) {
138     
139     cout<<"\n Event: "<<ev<<endl;
140     
141     // tracks from ITS
142     sprintf(trksName,"TreeT_ITS_%d",ev);
143     TTree *itsTree=(TTree*)itstrks->Get(trksName);
144     if(!itsTree) continue;
145     itsEntries = (Int_t)itsTree->GetEntries();
146     printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
147     
148     // tree from reference file
149     
150     sprintf(refsName,"Tree_Ref_%d",ev);
151     TTree *refTree=(TTree*)itsrefs->Get(refsName);
152     refTree->SetBranchAddress("rectracks",&rectrk);
153
154     //Getting primary Vertex
155     GetPrimaryVertex(ev,path);
156     
157     // count the total number of events
158     nTotEv++;
159
160     // call function which applies sigle track selection and
161     // separetes positives and negatives
162     TObjArray trksP(itsEntries/2);
163     Int_t *itsEntryP = new Int_t[itsEntries];
164     TObjArray trksN(itsEntries/2);
165     Int_t *itsEntryN = new Int_t[itsEntries];
166     SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN); 
167     
168     cout<<"#pos: "<<nTrksP<<endl;
169     cout<<"#neg: "<<nTrksN<<endl;
170     
171     // create the AliV0vertexer object
172     AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
173     
174     AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
175     
176     double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
177     
178     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
179       postrack = (AliITStrackV2*)trksP.At(iTrkP);
180
181       // get info on tracks PDG and mothers PDG from reference file
182       refTree->GetEvent(itsEntryP[iTrkP]);
183       pdg[0] = rectrk.pdg;
184       mum[0] = rectrk.mumpdg;
185       mumlab[0] = rectrk.mumlab;
186
187       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
188         negtrack = (AliITStrackV2*)trksN.At(iTrkN);
189
190         // get info on tracks PDG and mothers PDG from reference file
191         refTree->GetEvent(itsEntryN[iTrkN]);
192         pdg[1] = rectrk.pdg;
193         mum[1] = rectrk.mumpdg;
194         mumlab[1] = rectrk.mumlab;
195
196         D0->SetTracks(postrack,negtrack);
197
198         // ----------- DCA MINIMIZATION ------------------
199         //
200         // find the DCA and propagate the tracks to the DCA 
201         double dca = vertexer2->PropagateToDCA(negtrack,postrack);
202         
203         if(dca<cutDCA){
204           // define the AliV0vertex object
205           AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
206           // get position of the vertex
207           vertex2->GetXYZ(v2[0],v2[1],v2[2]);
208           delete vertex2;       
209           //if(D0->FindV0offline(v2) && D0->d0d0()){
210           if(D0->d0d0()){
211             D0->SetV0(v2);
212             D0->FindMomentaOffline();
213             if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
214               nD0++;
215               event[index]=ev; index++;
216               if(h){
217                 h1->Fill(D0->Pt());
218                 h3->Fill(D0->Eta());
219               }
220               
221               if(mumlab[0]==mumlab[1] && TMath::Abs(mum[0])==421 && TMath::Abs(mum[1])==421) { 
222                 nD0recSgn++;
223               } 
224               else { 
225                 if(TMath::Abs(mum[0])==421 || TMath::Abs(mum[0])==411 ||
226                    TMath::Abs(mum[1])==421 || TMath::Abs(mum[1])==411) {
227                   nD0recBkgS++;
228                 } else {
229                   nD0recBkg++;
230                 }  
231               }
232             }
233           } 
234         }
235       }
236     } 
237     
238     //delete D0;
239     //delete itstrks;
240     //delete itsTree;
241     //delete trksP;
242     //delete itsEntryP;
243     //delete trksN;
244     //delete itsEntryN;
245   }
246   
247   cout<<"\nMy #D0: "<<nD0<<"\n"<<endl;
248   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
249   printf("\n+++\n+++ Total number of D0 candidates:\n
250                  +++    Sgn:   %d\n
251                  +++    BkgS:  %d\n
252                  +++    Bkg:   %d\n+++\n",nD0recSgn,nD0recBkgS,nD0recBkg);
253
254   gBenchmark->Stop(name);
255   gBenchmark->Show(name);
256   
257   if(h){
258     if(!PtD0){
259       TCanvas *c = new TCanvas("c","c",700,1000);
260       c->Divide(1,2);
261       c->cd(1);
262       h1->Draw();
263       pt = new TPaveText(7.3,0.4,11,1.8,"br");
264       pt->SetFillColor(0);
265       pt->SetBorderSize(1);
266       pt->AddText("Cuts:");
267       Char_t st[1024];
268       sprintf(st,"First Pt:         %g",kPtCut);
269       pt->AddText(st);
270       sprintf(st,"d0 low:           %g",kd0Cut);
271       pt->AddText(st);
272       sprintf(st,"d0 high:          %g",kd0CutHigh);
273       pt->AddText(st);
274       sprintf(st,"V0 low:           %g",cuts[0]);
275       pt->AddText(st);
276       sprintf(st,"V0 high:          %g",cuts[1]);
277       pt->AddText(st);
278       sprintf(st,"InvMass Diff:     %g",cuts[2]);
279       pt->AddText(st);
280       sprintf(st,"cosPointingAngle: %g",cuts[3]);
281       pt->AddText(st);
282       sprintf(st,"d0d0:             %g",cuts[4]);
283       pt->AddText(st);
284       sprintf(st,"cosThetaStar:     %g",cuts[5]);
285       pt->AddText(st);
286       sprintf(st,"PtChild:          %g",cuts[6]);
287       pt->AddText(st);
288       sprintf(st,"DCA:              %g",cutDCA);
289       pt->AddText(st);
290       pt->Draw();
291       c_1->Modified();
292       c->cd();
293       c->cd(2);
294       h3->Draw();
295     }
296     if(PtD0){
297       PtD0(path,h2,h4);
298       TCanvas *c = new TCanvas("c","c",1000,700);
299       c->Divide(2,2);
300       c->cd(1);
301       h1->Draw();
302       pt = new TPaveText(7.3,0.4,11,1.8,"br");
303       pt->SetFillColor(0);
304       pt->SetBorderSize(1);
305       pt->AddText("Cuts:");
306       Char_t st[1024];
307       sprintf(st,"First Pt:         %g",kPtCut);
308       pt->AddText(st);
309       sprintf(st,"d0 low:           %g",kd0Cut);
310       pt->AddText(st);
311       sprintf(st,"d0 high:          %g",kd0CutHigh);
312       pt->AddText(st);
313       sprintf(st,"V0 low:           %g",cuts[0]);
314       pt->AddText(st);
315       sprintf(st,"V0 high:          %g",cuts[1]);
316       pt->AddText(st);
317       sprintf(st,"InvMass Diff:     %g",cuts[2]);
318       pt->AddText(st);
319       sprintf(st,"cosPointAng: %g",cuts[3]);
320       pt->AddText(st);
321       sprintf(st,"d0d0:             %g",cuts[4]);
322       pt->AddText(st);
323       sprintf(st,"cosTheta*:     %g",cuts[5]);
324       pt->AddText(st);
325       sprintf(st,"PtChild:          %g",cuts[6]);
326       pt->AddText(st);
327       sprintf(st,"DCA:              %g",cutDCA);
328       pt->AddText(st);
329       pt->Draw();
330       c_1->Modified();
331       c->cd();
332       c->cd(2);
333       h3->Draw();
334       c->cd(3);
335       h2->Draw();
336       c->cd(4);
337       h4->Draw();
338     }    
339   }
340   if(h){
341     if(!PtD0){
342       Char_t outName[1024];
343       sprintf(outName,"%s/ReconstructedD0.root",path);
344       TFile* outroot = new TFile(outName,"recreate");
345       h1->Write();
346       h3->Write();
347       outroot->Close();
348       delete outroot;
349     }
350   }
351   if(PtD0){
352     Char_t outName[1024];
353     sprintf(outName,"%s/ReconstructedD0.root",path);
354     TFile* outroot = new TFile(outName,"recreate");
355     h1->Write();
356     h2->Write();
357     h3->Write();
358     h4->Write();
359     outroot->Close();
360     delete outroot;
361   }
362   Char_t foutName[1024];
363   sprintf(foutName,"%s/Cuts",path);
364   ofstream fout(foutName);
365   Char_t st2[1024];
366   sprintf(st2,"First Pt:       %g",kPtCut);
367   fout<<st2<<endl;
368   sprintf(st2,"d0 low:         %g",kd0Cut);
369   fout<<st2<<endl;
370   sprintf(st2,"d0 high:        %g",kd0CutHigh);
371   fout<<st2<<endl;
372   sprintf(st2,"V0 low:         %g",cuts[0]);
373   fout<<st2<<endl;
374   sprintf(st2,"V0 high:        %g",cuts[1]);
375   fout<<st2<<endl;
376   sprintf(st2,"InvMass Diff:   %g",cuts[2]);
377   fout<<st2<<endl;
378   sprintf(st2,"cosPointAng:    %g",cuts[3]);
379   fout<<st2<<endl;
380   sprintf(st2,"d0d0:           %g",cuts[4]);
381   fout<<st2<<endl;
382   sprintf(st2,"cosTheta*:      %g",cuts[5]);
383   fout<<st2<<endl;
384   sprintf(st2,"PtChild:        %g",cuts[6]);
385   fout<<st2<<endl;
386   sprintf(st2,"DCA:            %g",cutDCA);
387   fout<<st2<<endl;
388   fout.close();
389
390   Char_t fName[1024];
391   sprintf(fName,"%s/Events",path);
392   ofstream fevent(fName);
393   for(int i=0;i<nD0;i++){fevent<<event[i]<<endl;}
394   fevent.close();
395 }
396 //___________________________________________________________________________
397 void   SelectTracks(TTree& itsTree,
398                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
399                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
400 //
401 // this function creates two TObjArrays with positive and negative tracks
402 //
403   nTrksP=0,nTrksN=0;
404
405  
406   Int_t entr = (Int_t)itsTree.GetEntries();
407
408   // trasfer tracks from tree to arrays
409   for(Int_t i=0; i<entr; i++) {
410
411     AliITStrackV2 *itstrack = new AliITStrackV2; 
412     itsTree.SetBranchAddress("tracks",&itstrack);
413
414     itsTree.GetEvent(i);
415
416     // single track selection
417     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
418
419     if(itstrack->Get1Pt()>0.) { // negative track
420       trksN.AddLast(itstrack);
421       itsEntryN[nTrksN] = i;
422       nTrksN++;
423     } else {                    // positive track
424       trksP.AddLast(itstrack);
425       itsEntryP[nTrksP] = i;
426       nTrksP++;
427     }
428
429   }
430
431   return;
432 }
433 //____________________________________________________________________________
434 Bool_t TrkCuts(const AliITStrackV2& trk) {
435 // 
436 // this function tells if track passes some kinematical cuts  
437 //
438   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
439   if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
440   //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
441
442   return kTRUE;
443 }
444 //____________________________________________________________________________
445 void GetPrimaryVertex(int i,Char_t* path="./") {
446
447   int event=i;
448
449   Char_t falice[1024];
450   sprintf(falice,"%s/galice.root",path);
451   TFile * galice = new TFile(falice);
452   
453   TDirectory * curdir;  
454
455   Char_t vname[20];
456   galice->cd();
457   
458   sprintf(vname,"Vertex_%d",event);
459   TArrayF o = 0;
460   o.Set(3);
461   AliHeader * header = 0;
462   
463   TTree * treeE = (TTree*)gDirectory->Get("TE");
464   treeE->SetBranchAddress("Header",&header);
465   treeE->GetEntry(event);
466   AliGenEventHeader* genHeader = header->GenEventHeader();
467   if(genHeader){
468     genHeader->PrimaryVertex(o);
469     primaryvertex[0] = (Double_t)o[0];
470     primaryvertex[1] = (Double_t)o[1];
471     primaryvertex[2] = (Double_t)o[2];
472   }
473   else{
474     printf("Can't find Header");
475   }
476   delete header;
477   delete galice;
478 }
479 //____________________________________________________________________________
480 PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
481   
482   Char_t falice[1024];
483   sprintf(falice,"%s/galice.root",path);  
484   TFile *f = new TFile(falice);
485   gAlice=(AliRun*)f->Get("gAlice");
486   
487   TParticle *p;
488   int nd0=0;
489   int nkminus =0;
490   int npipluss = 0;
491   int nEvent=gAlice->GetEventsPerRun();
492   
493   for (Int_t i = 0; i <nEvent; i++) {
494     cout<<"Event: "<<i<<endl;
495     gAlice->GetEvent(i);
496     Int_t nPart = gAlice->GetNtrack();
497     for (Int_t iPart = 0; iPart < nPart; iPart++) {
498       //cout<<"Particlenr.: "<<iPart<<endl;
499       p = (TParticle*)gAlice->Particle(iPart);
500       if (p->GetPdgCode()==421){
501         if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
502         h4->Fill(p.Eta());
503         nd0++;
504       }
505       if (p->GetPdgCode()==-321){
506         nkminus++;
507       }
508       if (p->GetPdgCode()==211){
509         npipluss++;
510       }
511     }
512   }
513   delete gAlice;
514 }