(martin) pt vs eta correction matrix calculation macro using CorrectionMatrix2D class.
[u/mrichter/AliRoot.git] / PWG0 / ptSpectra / makeCorrectionPtEta.C
CommitLineData
559d0c87 1//
2// Script to calculate PtvsEta correction map using the CorrectionMatrix2D class.
3//
4
5makeCorrectionPtEta(Char_t* dataDir, Int_t nRuns=20) {
6
7 Char_t str[256];
8
9 //gSystem->Load("/home/poghos/CERN2006/pp/esdTrackCuts/libESDtrackQuality.so");
10 gSystem->Load("libPWG0base");
11 gSystem->Load("libCorrectionMatrix2D.so");
12 // ########################################################
13 // selection of esd tracks
14 AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();
15 esdTrackCuts->DefineHistograms(1);
16 esdTrackCuts->SetMinNClustersTPC(50);
17 esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
18 esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
19 esdTrackCuts->SetRequireTPCRefit(kTRUE);
20 esdTrackCuts->SetMinNsigmaToVertex(3);
21 esdTrackCuts->SetAcceptKingDaughters(kFALSE);
22 AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
23 AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
24
25 // ########################################################
26 // definition of PtEta correction object
27 CorrectionMatrix2D* PtEtaMap = new CorrectionMatrix2D("CorrectionMatrix2");
28 PtEtaMap->SetHist("PtEta",80,0.,10.,120,-2.,2.);
29 //Float_t x[]={-20., -15., -5., 0., 6., 20.};
30 //Float_t y[]={-6. , -3., -1., 2., 3., 6. };
31 //PtEtaMap->SetHist("PtEta",5,x,5,y);
32 PtEtaMap->SetHistTitle("p_{t}","#eta");
33
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 // definition of used pointers
47 TFile* esdFile;
48 TTree* esdTree;
49 TBranch* esdBranch;
50
51 AliESD* esd =0;
52 // ########################################################
53 // loop over runs
54 Int_t nRunCounter = 0;
55 for (Int_t r=0; r<nDirs; r++) {
56
57 TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
58 if (!presentDir || !presentDir->IsDirectory())
59 continue;
60 // check that the files are there
61 TString currentDataDir;
62 currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
63 cout << "Processing directory " << currentDataDir.Data() << endl;
64 if ((!gSystem->Which(currentDataDir,"galice.root")) ||
65 (!gSystem->Which(currentDataDir,"AliESDs.root")))
66 continue;
67
68 if (nRunCounter++ >= nRuns)
69 break;
70
71 cout << "run #" << nRunCounter << endl;
72
73 // #########################################################
74 // setup galice and runloader
75 if (gAlice) {
76 delete gAlice->GetRunLoader();
77 delete gAlice;
78 gAlice=0;
79 }
80
81 sprintf(str,"%s/galice.root",currentDataDir.Data());
82 AliRunLoader* runLoader = AliRunLoader::Open(str);
83
84 runLoader->LoadgAlice();
85 gAlice = runLoader->GetAliRun();
86 runLoader->LoadKinematics();
87 runLoader->LoadHeader();
88
89
90 // #########################################################
91 // open esd file and get the tree
92
93 sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
94 // close it first to avoid memory leak
95 if (esdFile)
96 if (esdFile->IsOpen())
97 esdFile->Close();
98
99 esdFile = TFile::Open(str);
100 esdTree = (TTree*)esdFile->Get("esdTree");
101 if (!esdTree)
102 continue;
103 esdBranch = esdTree->GetBranch("ESD");
104 esdBranch->SetAddress(&esd);
105 if (!esdBranch)
106 continue;
107
108 // ########################################################
109 // Magnetic field
110 AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
111 //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
112
113 // ########################################################
114 // getting number of events
115 Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
116 Int_t nESDEvents = esdBranch->GetEntries();
117
118 if (nEvents!=nESDEvents) {
119 cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
120 return;
121 }
122 // ########################################################
123 // loop over number of events
124 cout << " looping over events..." << endl;
125 for(Int_t i=0; i<nEvents; i++) {
126
127 esdBranch->GetEntry(i);
128 runLoader->GetEvent(i);
129
130 // ########################################################
131 // get the EDS vertex
132 AliESDVertex* vtxESD = esd->GetVertex();
133
134 Double_t vtx[3];
135 Double_t vtx_res[3];
136 vtxESD->GetXYZ(vtx);
137
138 vtx_res[0] = vtxESD->GetXRes();
139 vtx_res[1] = vtxESD->GetYRes();
140 vtx_res[2] = vtxESD->GetZRes();
141
142 // the vertex should be reconstructed
143 if (strcmp(vtxESD->GetName(),"default")==0)
144 continue;
145
146 // the resolution should be reasonable???
147 if (vtx_res[2]==0 || vtx_res[2]>0.1)
148 continue;
149
150 // ########################################################
151 // loop over mc particles
152 AliStack* particleStack = runLoader->Stack();
153 Int_t nPrim = particleStack->GetNprimary();
154
155 for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
156
157 TParticle* part = particleStack->Particle(i_mc);
158 if (!part || strcmp(part->GetName(),"XXX")==0)
159 continue;
160
161 TParticlePDG* pdgPart = part->GetPDG();
162
163 Bool_t prim = kFALSE;
164 // check if it's a primary particle - is there a better way ???
165 if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
166 if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
167 prim = kTRUE;
168 }
169 if (!prim)
170 continue;
171
172 if (pdgPart->Charge()==0)
173 continue;
174
175 if (prim)
176 PtEtaMap->FillGene(part->Pt(), part->Eta());
177
178 }// end of mc particle
179
180 // ########################################################
181 // loop over esd tracks
182 Int_t nTracks = esd->GetNumberOfTracks();
183 for (Int_t t=0; t<nTracks; t++) {
184 AliESDtrack* esdTrack = esd->GetTrack(t);
185
186 // cut the esd track?
187 if (!esdTrackCuts->AcceptTrack(esdTrack))
188 continue;
189
190 Double_t Ptr0[3];
191 if(!esdTrack->GetPxPyPz(Ptr0)) continue;
192 fTrP= new TVector3(Ptr0);
193
194 PtEtaMap->FillMeas(fTrP->Pt(), fTrP->Eta());
195
196 } // end of track loop
197 } // end of event loop
198 } // end of run loop
199
200 PtEtaMap->Finish();
201
202 TFile* fout = new TFile("PtEtaCorrectionMap.root","RECREATE");
203
204 esdTrackCuts->SaveHistograms("EsdTrackCuts");
205 PtEtaMap->SaveHistograms();
206
207 fout->Write();
208 fout->Close();
209
210}