changed calculation of n_sigma to vertex (again...).
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / testESDtrackCuts.C
1 /* $Id$ */
2
3 testESDtrackCuts(Char_t* dataDir=, Int_t nRuns=10) {
4
5   Char_t str[256];
6
7   gSystem->Load("libESDtrackCuts.so");
8
9   // ########################################################
10   // definition of ESD track cuts
11
12   AliESDtrackCuts* trackCuts = new AliESDtrackCuts();
13   trackCuts->DefineHistograms(4);
14
15   trackCuts->SetMinNClustersTPC(50);
16   trackCuts->SetMaxChi2PerClusterTPC(3.5);
17   trackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
18   trackCuts->SetRequireTPCRefit(kTRUE);
19
20   trackCuts->SetMinNsigmaToVertex(3);
21   trackCuts->SetAcceptKingDaughters(kFALSE);
22
23   trackCuts->SetPRange(0.3);
24
25   AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
26
27   // ########################################################
28   // definition of used pointers
29   TFile* esdFile;
30   TTree* esdTree;
31   TBranch* esdBranch;
32
33   AliESD* esd = 0;
34
35   // ########################################################
36   // get the data dir  
37   Char_t execDir[256];
38   sprintf(execDir,gSystem->pwd());
39   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
40   TList* dirList            = baseDir->GetListOfFiles();
41   Int_t nDirs               = dirList->GetEntries();
42   // go back to the dir where this script is executed
43   gSystem->cd(execDir);
44   
45   // ########################################################
46   // loop over runs
47   Int_t nRunCounter = 0;
48   for (Int_t r=1; r<=nDirs; r++) {
49
50     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
51     if (!presentDir->IsDirectory())
52       continue;
53     // first check that the files are there
54     sprintf(str,"%s/%s",dataDir, presentDir->GetName());
55     if ((!gSystem->Which(str,"galice.root")) ||
56         (!gSystem->Which(str,"AliESDs.root"))) 
57       continue;
58     
59     if (nRunCounter++ >= nRuns)
60       break;    
61     
62     cout << "run #" << nRunCounter << endl;
63
64     // #########################################################
65     // setup galice and runloader
66     if (gAlice) {
67       delete gAlice->GetRunLoader();
68       delete gAlice;
69       gAlice=0;
70     }
71
72     sprintf(str,"%s/run%d/galice.root",dataDir,r);
73     AliRunLoader* runLoader = AliRunLoader::Open(str);
74
75     runLoader->LoadgAlice();
76     gAlice = runLoader->GetAliRun();
77     runLoader->LoadHeader();
78
79     // #########################################################
80     // open esd file and get the tree
81
82     sprintf(str,"%s/run%d/AliESDs.root",dataDir,r);
83     // close it first to avoid memory leak
84     if (esdFile)
85       if (esdFile->IsOpen())
86         esdFile->Close();
87
88     esdFile = TFile::Open(str);
89     esdTree = (TTree*)esdFile->Get("esdTree");
90     if (!esdTree)
91       continue;
92     esdBranch = esdTree->GetBranch("ESD");
93     esdBranch->SetAddress(&esd);
94     if (!esdBranch)
95       continue;
96
97     // ########################################################
98     // Magnetic field
99     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
100
101     // ########################################################
102     // getting number of events
103     Int_t nEvents    = (Int_t)runLoader->GetNumberOfEvents();
104     Int_t nESDEvents = esdBranch->GetEntries();
105     
106     if (nEvents!=nESDEvents) 
107       cout << " Warning: Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
108     
109     // ########################################################
110     // loop over number of events
111     cout << " looping over events..." << endl;
112     for(Int_t i=1; i<nEvents; i++) {
113       
114       esdBranch->GetEntry(i);
115       runLoader->GetEvent(i);
116       
117       // ########################################################
118       // get the EDS vertex
119       AliESDVertex* vtxESD = esd->GetVertex();
120       
121       Double_t vtxSigma[3];
122       vtxESD->GetSigmaXYZ(vtxSigma);
123       
124       // ########################################################
125       // loop over esd tracks      
126       Int_t nTracks = esd->GetNumberOfTracks();      
127
128       for (Int_t t=0; t<nTracks; t++) {
129         AliESDtrack* esdTrack = esd->GetTrack(t);      
130         
131         //trackCuts->AcceptTrack(esdTrack, vtxESD, esd->GetMagneticField());
132         trackCuts->AcceptTrack(esdTrack);
133         
134       } // end of track loop
135     } // end  of event loop
136   } // end of run loop
137
138   TFile* fout = new TFile("out.root","RECREATE");
139
140   trackCuts->SaveHistograms("esd_track_cuts");
141   
142   fout->Write();
143   fout->Close();
144
145 }