o) Added DrawHistograms function
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / testAnalysis.C
CommitLineData
75ec0f41 1//
2// Script to test the dN/dEta analysis using the dNdEtaAnalysis and
3// dNdEtaCorrection classes. Note that there is a cut on the events,
4// so the measurement will be biassed.
5//
6
7testAnalysis(Char_t* dataDir, Int_t nRuns=20) {
8
9 Char_t str[256];
10
11 gSystem->Load("../esdTrackCuts/libESDtrackQuality.so");
12 gSystem->Load("libdNdEta.so");
13
14 // ########################################################
15 // selection of esd tracks
16 ESDtrackQualityCuts* esdTrackCuts = new ESDtrackQualityCuts();
17 esdTrackCuts->DefineHistograms(1);
18
19 esdTrackCuts->SetMinNClustersTPC(50);
20 esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
21 esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
22 esdTrackCuts->SetRequireTPCRefit(kTRUE);
23
24 esdTrackCuts->SetMinNsigmaToVertex(3);
25 esdTrackCuts->SetAcceptKingDaughters(kFALSE);
26
27 AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
28
29 // ########################################################
30 // definition of dNdEta objects
31
32 dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection();
33
34 dNdEtaMap->LoadHistograms("correction_map.root","dndeta_correction");
35 dNdEtaMap->RemoveEdges(2,0,2);
36
37 dNdEtaAnalysis* analyse = new dNdEtaAnalysis("dndeta");
38
39 // ########################################################
40 // get the data dir
41 Char_t execDir[256];
42 sprintf(execDir,gSystem->pwd());
43 TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
44 TList* dirList = baseDir->GetListOfFiles();
45 Int_t nDirs = dirList->GetEntries();
46 // go back to the dir where this script is executed
47 gSystem->cd(execDir);
48
49
50 // ########################################################
51 // definition of used pointers
52 TFile* esdFile;
53 TTree* esdTree;
54 TBranch* esdBranch;
55
56 AliESD* esd =0;
57
58 // ########################################################
59 // loop over runs
60 Int_t nRunCounter = 0;
61 for (Int_t r=1; r<=nDirs; r++) {
62
63 TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
64 if (!presentDir->IsDirectory())
65 continue;
66 // check that the files are there
ceb5d1b5 67 TString currentDataDir;
68 currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
69 if ((!gSystem->Which(currentDataDir.Data(),"galice.root")) ||
70 (!gSystem->Which(currentDataDir.Data(),"AliESDs.root")))
75ec0f41 71 continue;
72
73 if (nRunCounter++ >= nRuns)
74 break;
75
76 cout << "run #" << nRunCounter << endl;
77
78 // #########################################################
79 // setup galice and runloader
80 if (gAlice) {
81 delete gAlice->GetRunLoader();
82 delete gAlice;
83 gAlice=0;
84 }
85
ceb5d1b5 86 sprintf(str,"%s/galice.root",currentDataDir.Data());
75ec0f41 87 AliRunLoader* runLoader = AliRunLoader::Open(str);
88
89 runLoader->LoadgAlice();
90 gAlice = runLoader->GetAliRun();
91 runLoader->LoadKinematics();
92 runLoader->LoadHeader();
93
94 // #########################################################
95 // open esd file and get the tree
96
ceb5d1b5 97 sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
75ec0f41 98 // close it first to avoid memory leak
99 if (esdFile)
100 if (esdFile->IsOpen())
101 esdFile->Close();
102
103 esdFile = TFile::Open(str);
104 esdTree = (TTree*)esdFile->Get("esdTree");
105 if (!esdTree)
106 continue;
107 esdBranch = esdTree->GetBranch("ESD");
108 esdBranch->SetAddress(&esd);
109 if (!esdBranch)
110 continue;
111
112 // ########################################################
113 // Magnetic field
114 AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
115 //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
116
117 // ########################################################
118 // getting number of events
119 Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
120 Int_t nESDEvents = esdBranch->GetEntries();
121
122 if (nEvents!=nESDEvents) {
123 cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
124 return;
125 }
126 // ########################################################
127 // loop over number of events
128 cout << " looping over events..." << endl;
129 for(Int_t i=1; i<nEvents; i++) {
130
131 esdBranch->GetEntry(i);
132 runLoader->GetEvent(i);
133
134
135 // ########################################################
136 // get the EDS vertex
137 AliESDVertex* vtxESD = esd->GetVertex();
138
139 Double_t vtx[3];
140 Double_t vtx_res[3];
141 vtxESD->GetXYZ(vtx);
142
143 vtx_res[0] = vtxESD->GetXRes();
144 vtx_res[1] = vtxESD->GetYRes();
145 vtx_res[2] = vtxESD->GetZRes();
146
147 // we need a good vertex
148 // => there will be a bias on the measurement, since this cuts away some events
149 if (strcmp(vtxESD->GetName(),"default")==0)
150 continue;
151 if (vtx_res[2]==0 || vtx_res[2]>0.1)
152 continue;
153
154 // ########################################################
155 // loop over esd tracks
156 Int_t nTracks = esd->GetNumberOfTracks();
157 for (Int_t t=0; t<nTracks; t++) {
158 AliESDtrack* esdTrack = esd->GetTrack(t);
159
160 if (!esdTrackCuts->AcceptTrack(esdTrack))
161 continue;
162
163 AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
164
165 if (tpcTrack->GetAlpha()==0) {
166 cout << " Warning esd track alpha = 0" << endl;
167 continue;
168 }
169
170 Float_t theta = tpcTrack->Theta();
171 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
172 Float_t correction = dNdEtaMap->GetCorrection(vtx[2], eta);
173
174 dNdEtaMap->FillMeas(vtx[2], eta);
175
176 analyse ->FillTrack(vtx[2], eta, correction);
177
178 } // end of track loop
179 analyse->FillEvent(vtx[2]);
180
181 } // end of event loop
182 } // end of run loop
183
184 analyse->Finish();
185
186 TFile* fout = new TFile("out.root","RECREATE");
187
188 esdTrackCuts->SaveHistograms("esd_tracks_cuts");
189 dNdEtaMap->SaveHistograms();
190 analyse ->SaveHistograms();
191
192 fout->Write();
193 fout->Close();
194
195}