559d0c87 |
1 | // |
2 | // Script to calculate PtvsEta correction map using the CorrectionMatrix2D class. |
3 | // |
4 | |
5 | makeCorrectionPtEta(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 | } |