]>
Commit | Line | Data |
---|---|---|
f72219cb | 1 | |
2 | /************************************************************************** | |
3 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
4 | * * | |
5 | * Author: The ALICE Off-line Project. * | |
6 | * Contributors are mentioned in the code where appropriate. * | |
7 | * * | |
8 | * Permission to use, copy, modify and distribute this software and its * | |
9 | * documentation strictly for non-commercial purposes is hereby granted * | |
10 | * without fee, provided that the above copyright notice appears in all * | |
11 | * copies and that both the copyright notice and this permission notice * | |
12 | * appear in the supporting documentation. The authors make no claims * | |
13 | * about the suitability of this software for any purpose. It is * | |
14 | * provided "as is" without express or implied warranty. * | |
15 | **************************************************************************/ | |
16 | ||
17 | /* | |
18 | ||
19 | ||
20 | Basic calibration and QA class for the TPC gain calibration based on tracks from BEAM events. | |
21 | ||
22 | ||
23 | Send comments etc. to: A.Kalweit@gsi.de, marian.ivanov@cern.ch | |
24 | */ | |
25 | ||
26 | ||
27 | #include "Riostream.h" | |
28 | #include "TChain.h" | |
29 | #include "TTree.h" | |
30 | #include "TH1F.h" | |
31 | #include "TH2F.h" | |
32 | #include "TList.h" | |
33 | #include "TMath.h" | |
34 | #include "TCanvas.h" | |
35 | #include "TFile.h" | |
36 | #include "TF1.h" | |
6feb400f | 37 | #include "TVectorF.h" |
f72219cb | 38 | #include "TProfile.h" |
39 | ||
40 | #include "AliTPCcalibDB.h" | |
41 | #include "AliTPCclusterMI.h" | |
42 | #include "AliTPCClusterParam.h" | |
43 | #include "AliTPCseed.h" | |
44 | #include "AliESDVertex.h" | |
45 | #include "AliESDEvent.h" | |
46 | #include "AliESDfriend.h" | |
47 | #include "AliESDInputHandler.h" | |
48 | #include "AliAnalysisManager.h" | |
49 | #include "AliTPCParam.h" | |
50 | ||
51 | #include "AliComplexCluster.h" | |
52 | #include "AliTPCclusterMI.h" | |
53 | ||
54 | #include "AliLog.h" | |
55 | ||
56 | #include "AliTPCcalibGainMult.h" | |
57 | ||
58 | #include "TTreeStream.h" | |
6feb400f | 59 | #include "TDatabasePDG.h" |
60 | #include "AliKFParticle.h" | |
61 | #include "AliKFVertex.h" | |
62 | #include "AliESDv0.h" | |
63 | #include "AliESDkink.h" | |
64 | #include "AliRecoParam.h" | |
65 | #include "AliTracker.h" | |
66 | #include "AliTPCTransform.h" | |
67 | #include "AliTPCROC.h" | |
68 | #include "TROOT.h" | |
f72219cb | 69 | |
70 | ClassImp(AliTPCcalibGainMult) | |
71 | ||
72 | ||
73 | AliTPCcalibGainMult::AliTPCcalibGainMult() | |
74 | :AliTPCcalibBase(), | |
75 | fMIP(0), | |
76 | fLowerTrunc(0), | |
77 | fUpperTrunc(0), | |
78 | fUseMax(kFALSE), | |
79 | fHistNTracks(0), | |
80 | fHistClusterShape(0), | |
81 | fHistQA(0), | |
82 | fHistGainSector(0), | |
83 | fHistPadEqual(0), | |
6feb400f | 84 | fHistGainMult(0), |
85 | fPIDMatrix(0), | |
86 | fHistdEdxMap(0), | |
87 | fHistdEdxMax(0), | |
88 | fHistdEdxTot(0), | |
89 | fdEdxTree(0), | |
90 | fBBParam(0) | |
f72219cb | 91 | { |
92 | // | |
93 | // Empty default cosntructor | |
94 | // | |
95 | AliInfo("Default Constructor"); | |
96 | } | |
97 | ||
98 | ||
99 | AliTPCcalibGainMult::AliTPCcalibGainMult(const Text_t *name, const Text_t *title) | |
100 | :AliTPCcalibBase(), | |
101 | fMIP(0), | |
102 | fLowerTrunc(0), | |
103 | fUpperTrunc(0), | |
104 | fUseMax(kFALSE), | |
105 | fHistNTracks(0), | |
106 | fHistClusterShape(0), | |
107 | fHistQA(0), | |
108 | fHistGainSector(0), | |
109 | fHistPadEqual(0), | |
6feb400f | 110 | fHistGainMult(0), |
111 | fPIDMatrix(0), | |
112 | fHistdEdxMap(0), | |
113 | fHistdEdxMax(0), | |
114 | fHistdEdxTot(0), | |
115 | fdEdxTree(0), | |
116 | fBBParam(0) | |
f72219cb | 117 | { |
118 | // | |
119 | // | |
120 | // | |
121 | SetName(name); | |
122 | SetTitle(title); | |
123 | // | |
124 | fMIP = 50.; | |
125 | fLowerTrunc = 0.02; // IMPORTANT CHANGE --> REMOVE HARDWIRED TRUNCATION FROM TRACKER | |
126 | fUpperTrunc = 0.6; | |
127 | fUseMax = kTRUE; // IMPORTANT CHANGE FOR PbPb; standard: kFALSE; | |
128 | // | |
129 | fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",1001,-0.5,1000.5); | |
130 | fHistClusterShape = new TH1F("fHistClusterShape","cluster shape; rms meas. / rms exp.;",300,0,3); | |
131 | fHistQA = new TH3F("fHistQA","dEdx; momentum p (GeV); TPC signal (a.u.); pad",500,0.1,20.,500,0.,500,6,-0.5,5.5); | |
6feb400f | 132 | AliTPCcalibBase::BinLogX(fHistQA); |
f72219cb | 133 | // |
134 | // | |
135 | // MIP, sect, pad (short,med,long,full,oroc), run, ncl | |
136 | Int_t binsGainSec[5] = { 145, 72, 4, 10000000, 65}; | |
137 | Double_t xminGainSec[5] = { 10., -0.5, -0.5, -0.5, -0.5}; | |
138 | Double_t xmaxGainSec[5] = {300., 71.5, 3.5, 9999999.5, 64.5}; | |
139 | TString axisNameSec[5]={"Q","sector","pad type","run", "ncl"}; | |
140 | TString axisTitleSec[5]={"Q (a.u)","sector","pad type","run","ncl"}; | |
141 | // | |
142 | fHistGainSector = new THnSparseF("fHistGainSector","0:MIP, 1:sect, 2:pad, 3:run, 4:ncl", 5, binsGainSec, xminGainSec, xmaxGainSec); | |
143 | for (Int_t iaxis=0; iaxis<5;iaxis++){ | |
144 | fHistGainSector->GetAxis(iaxis)->SetName(axisNameSec[iaxis]); | |
145 | fHistGainSector->GetAxis(iaxis)->SetTitle(axisTitleSec[iaxis]); | |
146 | } | |
147 | // | |
148 | // | |
149 | // | |
bf07be7e | 150 | Int_t binsPadEqual[6] = { 400, 400, 4, 20, 50, 100}; |
151 | Double_t xminPadEqual[6] = { 0.0, 0.0, -0.5, 0, -250, 0}; | |
152 | Double_t xmaxPadEqual[6] = { 2.0, 2.0, 3.5, 13000, 250, 3}; | |
f72219cb | 153 | TString axisNamePadEqual[6] = {"dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt"}; |
154 | TString axisTitlePadEqual[6] = {"dEdx_padRegion/mean_dEdx Qmax", "dEdx_padRegion/mean_dEdx Qtot","padType","mult","driftlength", "1/pt"}; | |
155 | // | |
156 | fHistPadEqual = new THnSparseF("fHistPadEqual","0:dEdx_pad/dEdx_mean, 1:pad, 2:mult, 3:drift, 4:1/pt", 6, binsPadEqual, xminPadEqual, xmaxPadEqual); | |
157 | for (Int_t iaxis=0; iaxis<6;iaxis++){ | |
158 | fHistPadEqual->GetAxis(iaxis)->SetName(axisNamePadEqual[iaxis]); | |
159 | fHistPadEqual->GetAxis(iaxis)->SetTitle(axisTitlePadEqual[iaxis]); | |
160 | } | |
161 | // | |
162 | // | |
163 | // MIP Qmax, MIP Qtot, z, pad, vtx. contribut., ncl | |
164 | Int_t binsGainMult[6] = { 145, 145, 25, 4, 100, 80}; | |
165 | Double_t xminGainMult[6] = { 10., 10., 0, -0.5, 0, -0.5}; | |
012c4694 | 166 | Double_t xmaxGainMult[6] = {300., 300., 250, 3.5, 13000, 159.5}; |
f72219cb | 167 | TString axisNameMult[6]={"Qmax","Qtot","drift","padtype""multiplicity","ncl"}; |
168 | TString axisTitleMult[6]={"Qmax (a.u)","Qtot (a.u.)","driftlenght l (cm)","Pad Type","multiplicity","ncl"}; | |
169 | // | |
170 | fHistGainMult = new THnSparseF("fHistGainMult","MIP Qmax, MIP Qtot, z, type, vtx. contribut.", 6, binsGainMult, xminGainMult, xmaxGainMult); | |
171 | for (Int_t iaxis=0; iaxis<6;iaxis++){ | |
172 | fHistGainMult->GetAxis(iaxis)->SetName(axisNameMult[iaxis]); | |
173 | fHistGainMult->GetAxis(iaxis)->SetTitle(axisTitleMult[iaxis]); | |
174 | } | |
175 | // | |
f72219cb | 176 | // |
6feb400f | 177 | // dedx maps - bigger granulatity in phi - |
178 | // to construct the dedx sector/phi map | |
179 | Int_t binsGainMap[4] = { 100, 90, 10, 6}; | |
180 | Double_t xminGainMap[4] = { 0.3, -TMath::Pi(), 0, 0}; | |
181 | Double_t xmaxGainMap[4] = { 2, TMath::Pi(), 1, 6}; | |
182 | TString axisNameMap[4] = {"Q_Qexp","phi", "1/Qexp","Pad Type"}; | |
183 | TString axisTitleMap[4] = {"Q/Q_{exp}","#phi (a.u.)","1/Q_{exp}","Pad Type"}; | |
184 | // | |
185 | fHistdEdxMap = new THnSparseF("fHistdEdxMap","fHistdEdxMap", 4, binsGainMap, xminGainMap, xmaxGainMap); | |
186 | for (Int_t iaxis=0; iaxis<4;iaxis++){ | |
187 | fHistdEdxMap->GetAxis(iaxis)->SetName(axisNameMap[iaxis]); | |
188 | fHistdEdxMap->GetAxis(iaxis)->SetTitle(axisTitleMap[iaxis]); | |
189 | } | |
190 | // | |
191 | // | |
192 | // | |
193 | // dedx maps | |
194 | Int_t binsGainMax[6] = { 100, 10, 10, 10, 5, 3}; | |
195 | Double_t xminGainMax[6] = { 0.5, 0, 0, 0, 0, 0}; | |
196 | Double_t xmaxGainMax[6] = { 1.5, 1, 1.0, 1.0, 3000, 3}; | |
197 | TString axisNameMax[6] = {"Q_Qexp","1/Qexp", "phi","theta","mult", "Pad Type"}; | |
198 | TString axisTitleMax[6] = {"Q/Q_{exp}","1/Qexp", "#phi","#theta","mult","Pad Type"}; | |
199 | // | |
200 | fHistdEdxMax = new THnSparseF("fHistdEdxMax","fHistdEdxMax", 6, binsGainMax, xminGainMax, xmaxGainMax); | |
201 | fHistdEdxTot = new THnSparseF("fHistdEdxTot","fHistdEdxTot", 6, binsGainMax, xminGainMax, xmaxGainMax); | |
202 | for (Int_t iaxis=0; iaxis<6;iaxis++){ | |
203 | fHistdEdxMax->GetAxis(iaxis)->SetName(axisNameMax[iaxis]); | |
204 | fHistdEdxMax->GetAxis(iaxis)->SetTitle(axisTitleMax[iaxis]); | |
205 | fHistdEdxTot->GetAxis(iaxis)->SetName(axisNameMax[iaxis]); | |
206 | fHistdEdxTot->GetAxis(iaxis)->SetTitle(axisTitleMax[iaxis]); | |
207 | } | |
208 | // | |
209 | AliInfo("Non Default Constructor"); | |
f72219cb | 210 | } |
211 | ||
212 | ||
213 | AliTPCcalibGainMult::~AliTPCcalibGainMult(){ | |
214 | // | |
6feb400f | 215 | // Destructor |
f72219cb | 216 | // |
b9ab8e40 | 217 | delete fHistNTracks; // histogram showing number of ESD tracks per event |
218 | delete fHistClusterShape; // histogram to check the cluster shape | |
219 | delete fHistQA; // dE/dx histogram showing the final spectrum | |
220 | // | |
221 | delete fHistGainSector; // histogram which shows MIP peak for each of the 3x36 sectors (pad region) | |
222 | delete fHistPadEqual; // histogram for the equalization of the gain in the different pad regions -> pass0 | |
223 | delete fHistGainMult; // histogram which shows decrease of MIP signal as a function | |
6feb400f | 224 | // |
225 | delete fHistdEdxMap; | |
226 | delete fHistdEdxMax; | |
227 | delete fHistdEdxTot; | |
228 | delete fdEdxTree; | |
f72219cb | 229 | } |
230 | ||
231 | ||
232 | ||
233 | void AliTPCcalibGainMult::Process(AliESDEvent *event) { | |
234 | // | |
6feb400f | 235 | // Main function of the class |
236 | // 1. Select Identified particles - for identified particles the flag in the PID matrix is stored | |
237 | // 1.a) ProcessV0s - select Electron (gamma coversion) and pion canditates (from K0s) | |
238 | // 1.b) ProcessTOF - select - Proton, kaon and pions candidates | |
239 | // AS THE TOF not calibrated yet in Pass0 - we are calibrating the TOF T0 in this function | |
240 | // 1.c) ProcessCosmic - select cosmic mumn candidates - too few entries - not significant for the calibration | |
241 | // 1.d) ProcessKinks - select Kaon and pion candidates. From our experience (see Kink debug streamer), the angular cut for kink daughter is not sufficient - big contamination - delta rays, hadronic interaction (proton) | |
242 | // - NOT USED for the | |
243 | // | |
244 | // 2. Loop over tracks | |
245 | // 2.a DumpTrack() - for identified particles dump the track and dEdx information into the tree (for later fitting) | |
246 | // 3. Actual fitting for the moment macro | |
247 | ||
f72219cb | 248 | // |
6feb400f | 249 | // Criteria for the track selection |
f72219cb | 250 | // |
6feb400f | 251 | const Int_t kMinNCL=80; // minimal number of cluster - tracks accepted for the dedx calibration |
252 | const Double_t kMaxEta=0.8; // maximal eta fo the track to be accepted | |
253 | const Double_t kMaxDCAR=10; // maximal DCA R of the track | |
254 | const Double_t kMaxDCAZ=5; // maximal DCA Z of the track | |
255 | const Double_t kMIPPt=0.45; // MIP pt | |
256 | ||
f72219cb | 257 | if (!event) { |
258 | Printf("ERROR: ESD not available"); | |
259 | return; | |
260 | } | |
6feb400f | 261 | fCurrentEvent=event ; |
262 | fMagF = event->GetMagneticField(); | |
263 | Int_t ntracks=event->GetNumberOfTracks(); | |
f72219cb | 264 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend")); |
265 | if (!esdFriend) { | |
6feb400f | 266 | //Printf("ERROR: esdFriend not available"); |
267 | delete fPIDMatrix; | |
268 | return; | |
269 | } | |
270 | if (!(esdFriend->TestSkipBit())) fPIDMatrix= new TMatrixD(ntracks,5); | |
271 | fHistNTracks->Fill(ntracks); | |
272 | ProcessCosmic(event); // usually not enogh statistic | |
273 | ||
274 | if (esdFriend->TestSkipBit()) { | |
275 | return; | |
f72219cb | 276 | } |
6feb400f | 277 | // |
278 | ProcessV0s(event); // | |
279 | ProcessTOF(event); // | |
280 | ProcessKinks(event); // not relyable | |
281 | DumpHPT(event); // | |
f72219cb | 282 | UInt_t runNumber = event->GetRunNumber(); |
012c4694 | 283 | Int_t nContributors = event->GetNumberOfTracks(); |
f72219cb | 284 | // |
285 | // track loop | |
286 | // | |
287 | for (Int_t i=0;i<ntracks;++i) { | |
288 | // | |
289 | AliESDtrack *track = event->GetTrack(i); | |
290 | if (!track) continue; | |
291 | // | |
292 | AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam(); | |
293 | if (!trackIn) continue; | |
294 | ||
295 | // calculate necessary track parameters | |
296 | Double_t meanP = trackIn->GetP(); | |
297 | Int_t ncls = track->GetTPCNcls(); | |
298 | ||
6feb400f | 299 | if (ncls < kMinNCL) continue; |
f72219cb | 300 | // exclude tracks which do not look like primaries or are simply too short or on wrong sectors |
6feb400f | 301 | if (TMath::Abs(trackIn->Eta()) > kMaxEta) continue; |
f72219cb | 302 | |
f72219cb | 303 | UInt_t status = track->GetStatus(); |
304 | if ((status&AliESDtrack::kTPCrefit)==0) continue; | |
305 | //if (track->GetNcls(0) < 3) continue; // ITS clusters | |
306 | Float_t dca[2], cov[3]; | |
307 | track->GetImpactParameters(dca,cov); | |
308 | Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0]); | |
6feb400f | 309 | if (primVtxDCA > kMaxDCAR || primVtxDCA < 0.00001) continue; |
310 | if (TMath::Abs(dca[1]) > kMaxDCAZ) continue; | |
f72219cb | 311 | // |
f72219cb | 312 | // |
6feb400f | 313 | // fill Alexander QA histogram |
f72219cb | 314 | // |
f72219cb | 315 | if (primVtxDCA < 3 && track->GetNcls(0) > 3 && track->GetKinkIndex(0) == 0 && ncls > 100) fHistQA->Fill(meanP, track->GetTPCsignal(), 5); |
316 | ||
317 | // Get seeds | |
318 | AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i); | |
319 | if (!friendTrack) continue; | |
320 | TObject *calibObject; | |
321 | AliTPCseed *seed = 0; | |
322 | for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { | |
323 | if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
324 | } | |
6feb400f | 325 | if (seed) DumpTrack(track, friendTrack,seed,i); // MI implementation for the identified particles |
326 | // | |
327 | if (seed) { // seed the container with track parameters and the clusters | |
328 | // | |
329 | const AliExternalTrackParam * trackOut = friendTrack->GetTPCOut(); // tack at the outer radius of TPC | |
f72219cb | 330 | if (!trackIn) continue; |
331 | if (!trackOut) continue; | |
332 | Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ()); | |
333 | // | |
334 | for (Int_t irow =0; irow<160;irow++) { | |
335 | AliTPCTrackerPoint * point = seed->GetTrackPoint(irow); | |
336 | if (point==0) continue; | |
337 | AliTPCclusterMI * cl = seed->GetClusterPointer(irow); | |
338 | if (cl==0) continue; | |
339 | // | |
340 | Float_t rsigmay = TMath::Sqrt(point->GetSigmaY()); | |
341 | fHistClusterShape->Fill(rsigmay); | |
342 | } | |
343 | // | |
344 | Int_t row0 = 0; | |
345 | Int_t row1 = 160; | |
346 | // | |
347 | Double_t signalShortMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,0,62); | |
348 | Double_t signalMedMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,63,126); | |
349 | Double_t signalLongMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,127,159); | |
350 | Double_t signalMax = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1); | |
351 | Double_t signalArrayMax[4] = {signalShortMax, signalMedMax, signalLongMax, signalMax}; | |
352 | // | |
353 | Double_t signalShortTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,0,62); | |
354 | Double_t signalMedTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,63,126); | |
355 | Double_t signalLongTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,127,159); | |
356 | Double_t signalTot = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1); | |
357 | Double_t signalArrayTot[4] = {signalShortTot, signalMedTot, signalLongTot, signalTot}; | |
358 | // | |
359 | Double_t mipSignalShort = fUseMax ? signalShortMax : signalShortTot; | |
360 | Double_t mipSignalMed = fUseMax ? signalMedMax : signalMedTot; | |
361 | Double_t mipSignalLong = fUseMax ? signalLongMax : signalLongTot; | |
362 | Double_t mipSignalOroc = seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax,63,159); | |
363 | Double_t signal = fUseMax ? signalMax : signalTot; | |
364 | // | |
365 | fHistQA->Fill(meanP, mipSignalShort, 0); | |
366 | fHistQA->Fill(meanP, mipSignalMed, 1); | |
367 | fHistQA->Fill(meanP, mipSignalLong, 2); | |
368 | fHistQA->Fill(meanP, signal, 3); | |
369 | fHistQA->Fill(meanP, mipSignalOroc, 4); | |
370 | // | |
371 | // "dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt" | |
372 | Float_t meanMax = (1/3.)*(signalArrayMax[0] + signalArrayMax[1] + signalArrayMax[2]); | |
b4ffdb90 | 373 | Float_t meanTot = (1/3.)*(signalArrayTot[0] + signalArrayTot[1] + signalArrayTot[2]); |
4a4c6be2 | 374 | if (meanMax < 1e-5 || meanTot < 1e-5) continue; |
f72219cb | 375 | for(Int_t ipad = 0; ipad < 4; ipad ++) { |
376 | Double_t vecPadEqual[6] = {signalArrayMax[ipad]/meanMax, signalArrayTot[ipad]/meanTot, ipad, nContributors, meanDrift, track->OneOverPt()}; | |
377 | fHistPadEqual->Fill(vecPadEqual); | |
378 | } | |
379 | // | |
6feb400f | 380 | // if (meanP > 0.4 && meanP < 0.55) { |
381 | if ( TMath::Abs(meanP-kMIPPt)<0.05 ) { | |
f72219cb | 382 | Double_t vecMult[6] = {seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,1,row0,row1), |
383 | seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,0,row0,row1), | |
384 | meanDrift, | |
385 | 3, | |
386 | nContributors, | |
387 | ncls}; | |
388 | // | |
389 | fHistGainMult->Fill(vecMult); | |
390 | vecMult[0]=mipSignalShort; vecMult[1]=mipSignalShort; vecMult[3]=0; | |
391 | fHistGainMult->Fill(vecMult); | |
392 | vecMult[0]=mipSignalMed; vecMult[1]=mipSignalMed; vecMult[3]=1; | |
393 | fHistGainMult->Fill(vecMult); | |
394 | vecMult[0]=mipSignalLong; vecMult[1]=mipSignalLong; vecMult[3]=2; | |
395 | fHistGainMult->Fill(vecMult); | |
396 | // | |
397 | } | |
398 | // | |
399 | // | |
6feb400f | 400 | if ( TMath::Abs(meanP-kMIPPt)>0.05 ) continue; // only MIP pions |
401 | //if (meanP > 0.5 || meanP < 0.4) continue; // only MIP pions | |
f72219cb | 402 | // |
403 | // for each track, we look at the three different pad regions, split it into tracklets, check if the sector does not change and fill the histogram | |
404 | // | |
405 | Bool_t isNotSplit[3] = {kTRUE, kTRUE, kTRUE}; // short, medium, long (true if the track is not split between two chambers) | |
406 | // | |
407 | Double_t sector[4] = {-1, -1, -1, -1}; // sector number short, medium, long, all | |
408 | Int_t ncl[3] = {0,0,0}; | |
409 | // | |
6feb400f | 410 | |
f72219cb | 411 | for (Int_t irow=0; irow < 159; irow++){ |
412 | Int_t padRegion = 0; | |
413 | if (irow > 62) padRegion = 1; | |
414 | if (irow > 126) padRegion = 2; | |
415 | // | |
416 | AliTPCclusterMI* cluster = seed->GetClusterPointer(irow); | |
417 | if (!cluster) continue; | |
418 | if (sector[padRegion] == -1) { | |
419 | sector[padRegion] = cluster->GetDetector(); | |
420 | continue; | |
421 | } | |
422 | if (sector[padRegion] != -1 && sector[padRegion] != cluster->GetDetector()) isNotSplit[padRegion] = kFALSE; | |
423 | ncl[padRegion]++; | |
424 | } | |
425 | // | |
426 | // MIP, sect, pad, run | |
427 | // | |
428 | Double_t vecMip[5] = {mipSignalShort, mipSignalMed, mipSignalLong, signal, mipSignalOroc}; | |
429 | // | |
430 | for(Int_t ipad = 0; ipad < 3; ipad++) { | |
6feb400f | 431 | // AK. - run Number To be removed - not needed |
f72219cb | 432 | Double_t vecGainSec[5] = {vecMip[ipad], sector[ipad], ipad, runNumber, ncl[ipad]}; |
433 | if (isNotSplit[ipad]) fHistGainSector->Fill(vecGainSec); | |
434 | } | |
435 | } | |
436 | ||
437 | } | |
6feb400f | 438 | |
439 | delete fPIDMatrix; | |
f72219cb | 440 | } |
441 | ||
442 | ||
a8ef8a9c | 443 | void AliTPCcalibGainMult::MakeLookup(THnSparse * /*hist*/, Char_t * /*outputFile*/) { |
444 | // | |
445 | // Not yet implemented | |
446 | // | |
4a4c6be2 | 447 | } |
f72219cb | 448 | |
449 | ||
450 | void AliTPCcalibGainMult::Analyze() { | |
6feb400f | 451 | // |
452 | // Not implemented | |
453 | // | |
f72219cb | 454 | |
455 | return; | |
456 | ||
457 | } | |
458 | ||
459 | ||
460 | Long64_t AliTPCcalibGainMult::Merge(TCollection *li) { | |
6feb400f | 461 | // |
462 | // merging of the component | |
463 | // | |
f72219cb | 464 | |
465 | TIterator* iter = li->MakeIterator(); | |
466 | AliTPCcalibGainMult* cal = 0; | |
467 | ||
468 | while ((cal = (AliTPCcalibGainMult*)iter->Next())) { | |
469 | if (!cal->InheritsFrom(AliTPCcalibGainMult::Class())) { | |
470 | Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName()); | |
471 | return -1; | |
472 | } | |
473 | ||
474 | if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks()); | |
475 | if (cal->GetHistClusterShape()) fHistClusterShape->Add(cal->GetHistClusterShape()); | |
476 | if (cal->GetHistQA()) fHistQA->Add(cal->GetHistQA()); | |
477 | if (cal->GetHistGainSector()) fHistGainSector->Add(cal->GetHistGainSector()); | |
478 | if (cal->GetHistPadEqual()) fHistPadEqual->Add(cal->GetHistPadEqual()); | |
479 | if (cal->GetHistGainMult()) fHistGainMult->Add(cal->GetHistGainMult()); | |
6feb400f | 480 | if (cal->fHistdEdxMap){ |
481 | if (fHistdEdxMap) fHistdEdxMap->Add(cal->fHistdEdxMap); | |
482 | } | |
483 | if (cal->fHistdEdxMax){ | |
484 | if (fHistdEdxMax) fHistdEdxMax->Add(cal->fHistdEdxMax); | |
485 | } | |
486 | if (cal->fHistdEdxTot){ | |
487 | if (fHistdEdxTot) fHistdEdxTot->Add(cal->fHistdEdxTot); | |
488 | } | |
489 | // | |
490 | // Originally we tireied to write the tree to the same file as other calibration components | |
491 | // We failed in merging => therefore this optio was disabled | |
492 | // | |
493 | // if (cal->fdEdxTree && cal->fdEdxTree->GetEntries()>0) { | |
494 | // if (fdEdxTree) { | |
495 | // const Int_t kMax=100000; | |
496 | // Int_t entriesSum = (Int_t)fdEdxTree->GetEntries(); | |
497 | // Int_t entriesCurrent = (Int_t)cal->fdEdxTree->GetEntries(); | |
498 | // Int_t entriesCp=TMath::Min((Int_t) entriesCurrent*(kMax*entriesSum),entriesCurrent); | |
499 | // // cal->fdEdxTree->SetBranchStatus("track*",0); | |
500 | // // cal->fdEdxTree->SetBranchStatus("vertex*",0); | |
501 | // // cal->fdEdxTree->SetBranchStatus("tpcOut*",0); | |
502 | // // cal->fdEdxTree->SetBranchStatus("vec*",0); | |
503 | // // fdEdxTree->SetBranchStatus("track*",0); | |
504 | // // fdEdxTree->SetBranchStatus("vertex*",0); | |
505 | // // fdEdxTree->SetBranchStatus("tpcOut*",0); | |
506 | // // fdEdxTree->SetBranchStatus("vec*",0); | |
507 | // fdEdxTree->Print(); | |
508 | // fdEdxTree->Dump(); | |
509 | // fdEdxTree->GetEntry(0); | |
510 | // for (Int_t i=0; i<entriesCurrent; i++){ | |
511 | // cal->fdEdxTree->CopyAddresses(fdEdxTree); | |
512 | // cal->fdEdxTree->GetEntry(i); | |
513 | // fdEdxTree->Fill(); | |
514 | // } | |
515 | // TObjArray *brarray = cal->fdEdxTree->GetListOfBranches(); | |
516 | // for (Int_t i=0; i<brarray->GetEntries(); i++) {TBranch * br = (TBranch *)brarray->At(i); br->SetAddress(0); } | |
517 | // } | |
518 | // else{ | |
519 | // fdEdxTree = (TTree*)(cal->fdEdxTree->Clone()); | |
520 | // TObjArray *brarray = fdEdxTree->GetListOfBranches(); | |
521 | // for (Int_t i=0; i<brarray->GetEntries(); i++) {TBranch * br = (TBranch *)brarray->At(i); br->SetAddress(0);} | |
522 | // } | |
523 | //} | |
524 | ||
f72219cb | 525 | } |
526 | ||
527 | return 0; | |
528 | ||
529 | } | |
530 | ||
531 | ||
532 | ||
f72219cb | 533 | |
534 | ||
535 | void AliTPCcalibGainMult::UpdateGainMap() { | |
536 | // | |
537 | // read in the old gain map and scale it appropriately... | |
538 | // | |
539 | /* | |
540 | gSystem->Load("libANALYSIS"); | |
541 | gSystem->Load("libTPCcalib"); | |
542 | // | |
543 | TFile jj("Run0_999999999_v1_s0.root"); | |
544 | AliTPCCalPad * pad = AliCDBEntry->GetObject()->Clone(); | |
545 | TFile hh("output.root"); | |
546 | AliTPCcalibGainMult * gain = calibTracksGain; | |
547 | TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,1); | |
548 | // | |
549 | TObjArray arr; | |
550 | histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr); | |
551 | TH1D * meanGainSec = arr->At(1); | |
552 | Double_t gainsIROC[36]; | |
553 | Double_t gainsOROC[36]; | |
554 | Double_t gains[72]; | |
555 | // | |
556 | for(Int_t isec = 1; isec < meanGainSec->GetNbinsX() + 1; isec++) { | |
557 | cout << isec << " " << meanGainSec->GetXaxis()->GetBinCenter(isec) << " " <<meanGainSec->GetBinContent(isec) << endl; | |
558 | gains[isec-1] = meanGainSec->GetBinContent(isec); | |
559 | if (isec < 37) { | |
560 | gainsIROC[isec-1] = meanGainSec->GetBinContent(isec); | |
561 | } else { | |
562 | gainsOROC[isec - 36 -1] = meanGainSec->GetBinContent(isec); | |
563 | } | |
564 | } | |
565 | Double_t meanIroc = TMath::Mean(36, gainsIROC); | |
566 | Double_t meanOroc = TMath::Mean(36, gainsIROC); | |
567 | for(Int_t i = 0; i < 36; i++) gains[i] /= meanIroc; | |
568 | for(Int_t i = 36; i < 72; i++) gains[i] /= meanOroc; | |
569 | // | |
570 | for(Int_t i = 0; i < 72; i++) { | |
571 | AliTPCCalROC * chamber = pad->GetCalROC(i); | |
572 | chamber->Multiply(gains[i]); | |
573 | cout << i << " "<< chamber->GetMean() << endl; | |
574 | } | |
575 | // | |
576 | // update the OCDB | |
577 | // | |
578 | AliCDBMetaData *metaData= new AliCDBMetaData(); | |
579 | metaData->SetObjectClassName("AliTPCCalPad"); | |
580 | metaData->SetResponsible("Alexander Kalweit"); | |
581 | metaData->SetBeamPeriod(1); | |
582 | metaData->SetAliRootVersion("04-19-05"); //root version | |
583 | metaData->SetComment("New gain map for 1600V OROC gain increase and equalization. Valid for runs starting after technical stop beginning of September."); | |
584 | AliCDBId id1("TPC/Calib/GainFactorDedx", 131541, AliCDBRunRange::Infinity()); // important: new gain runs here.. | |
585 | AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///d/alice05/akalweit/projects/OCDBupdate/HighGain_2010-09-03/OCDB/"); | |
586 | gStorage->Put(pad, id1, metaData); | |
587 | */ | |
588 | ||
589 | } | |
590 | ||
591 | void AliTPCcalibGainMult::UpdateClusterParam() { | |
592 | // | |
593 | // | |
594 | // | |
595 | /* | |
596 | gSystem->Load("libANALYSIS"); | |
597 | gSystem->Load("libTPCcalib"); | |
598 | // | |
599 | TFile ff("OldClsParam.root"); | |
600 | AliTPCClusterParam * param = AliCDBEntry->GetObject()->Clone(); | |
601 | ||
602 | TFile hh("output.root"); | |
603 | AliTPCcalibGainMult * gain = calibGainMult; | |
604 | TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,2); | |
605 | TObjArray arr; | |
606 | histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr); | |
607 | histGainSec->Draw("colz"); | |
608 | TH1D * fitVal = arr.At(1); | |
609 | fitVal->Draw("same"); | |
610 | param->GetQnormCorrMatrix()->Print(); | |
611 | param->GetQnormCorrMatrix()(0,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qtot | |
612 | param->GetQnormCorrMatrix()(1,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qtot | |
613 | param->GetQnormCorrMatrix()(2,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qtot | |
614 | // | |
615 | param->GetQnormCorrMatrix()(3,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qmax -> scaling assumed | |
616 | param->GetQnormCorrMatrix()(4,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qmax -> scaling assumed | |
617 | param->GetQnormCorrMatrix()(5,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qmax -> scaling assumed | |
618 | // | |
619 | TFile jj("newClusterParam.root","RECREATE"); | |
620 | param->Write(); | |
621 | param->GetQnormCorrMatrix()->Print(); | |
622 | // | |
623 | // update the OCDB | |
624 | // | |
625 | AliCDBMetaData *metaData= new AliCDBMetaData(); | |
626 | metaData->SetObjectClassName("AliTPCClusterParam"); | |
627 | metaData->SetResponsible("Alexander Kalweit"); | |
628 | metaData->SetBeamPeriod(1); | |
629 | metaData->SetAliRootVersion("04-19-04"); //root version | |
630 | metaData->SetComment("1600V OROC / hard thres. / new algorithm"); | |
631 | AliCDBId id1("TPC/Calib/ClusterParam", 0, AliCDBRunRange::Infinity()); // important: new gain runs here.. | |
632 | AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///lustre/alice/akalweit/baseline/CalibrationEntries/OldThres_NewAlgo_PP"); | |
633 | gStorage->Put(param, id1, metaData); | |
634 | */ | |
635 | ||
636 | ||
637 | } | |
638 | ||
6feb400f | 639 | |
640 | void AliTPCcalibGainMult::DumpTrack(AliESDtrack * track, AliESDfriendTrack *ftrack, AliTPCseed * seed, Int_t index){ | |
641 | // | |
642 | // dump interesting tracks | |
643 | // 1. track at MIP region | |
644 | // 2. highly ionizing protons | |
645 | // pidType: 0 - unselected | |
646 | // 1 - TOF | |
647 | // 2 - V0 | |
648 | // 4 - Cosmic | |
649 | // or of value | |
650 | // | |
651 | const Int_t kMax=10000; | |
652 | const Int_t kMinRows=80; | |
653 | const Double_t kDCAcut=30; | |
654 | // | |
655 | // Bethe Bloch paramerization | |
656 | // | |
657 | Double_t kp1= 0.0851148; | |
658 | Double_t kp2= 9.25771; | |
659 | Double_t kp3= 2.6558e-05; | |
660 | Double_t kp4= 2.32742; | |
661 | Double_t kp5= 1.83039; | |
662 | if (fBBParam){ | |
663 | kp1=(*fBBParam)[0]; | |
664 | kp2=(*fBBParam)[1]; | |
665 | kp3=(*fBBParam)[2]; | |
666 | kp4=(*fBBParam)[3]; | |
667 | kp5=(*fBBParam)[4]; | |
668 | } | |
669 | // | |
670 | AliTPCROC *roc = AliTPCROC::Instance(); | |
671 | TDatabasePDG *pdg = TDatabasePDG::Instance(); | |
672 | ||
673 | Int_t nclITS = track->GetNcls(0); | |
674 | Int_t ncl = track->GetTPCncls(); | |
675 | Double_t ncl21 = track->GetTPCClusterInfo(3,1); | |
676 | Double_t ncl20 = track->GetTPCClusterInfo(3,0); | |
677 | // | |
678 | if (!seed) return; | |
679 | if (ncl21<kMinRows) return; | |
680 | static Int_t counter=0; | |
681 | static Int_t counterHPT=0; | |
682 | // | |
683 | static TH1F *hisBB=new TH1F("hisBB","hisBB",20,0.1,1.00); // bethe bloch histogram = | |
684 | // used to cover more homegenously differnt dEdx regions | |
685 | static Double_t massPi = pdg->GetParticle("pi-")->Mass(); // | |
686 | static Double_t massK = pdg->GetParticle("K-")->Mass(); | |
687 | static Double_t massP = pdg->GetParticle("proton")->Mass(); | |
688 | static Double_t massE = pdg->GetParticle("e-")->Mass(); | |
689 | static Double_t massMuon = pdg->GetParticle("mu-")->Mass(); | |
690 | static Double_t radius0= roc->GetPadRowRadiiLow(roc->GetNRows(0)/2); | |
691 | static Double_t radius1= roc->GetPadRowRadiiUp(30); | |
692 | static Double_t radius2= roc->GetPadRowRadiiUp(roc->GetNRows(36)-15); | |
693 | ||
694 | AliESDVertex *vertex= (AliESDVertex *)fCurrentEvent->GetPrimaryVertex(); | |
695 | // | |
696 | // Estimate current MIP position - | |
697 | // | |
698 | static Double_t mipArray[kMax]; // mip array | |
699 | static Int_t counterMIP0=0; | |
700 | static Double_t medianMIP0=100000; // current MIP median position - estimated after some mimnimum number of MIP tracks | |
701 | ||
702 | if (TMath::Abs(track->GetP()-0.5)<0.1&&track->GetTPCsignal()/medianMIP0<1.5){ | |
703 | mipArray[counterMIP0%kMax]= track->GetTPCsignal(); | |
704 | counterMIP0++; | |
705 | if (counterMIP0>10) medianMIP0=TMath::Median(counterMIP0%kMax, mipArray); | |
706 | } | |
707 | // the PID as defiend from the external sources | |
708 | // | |
709 | Int_t isElectron = TMath::Nint((*fPIDMatrix)(index,0)); | |
710 | Int_t isMuon = TMath::Nint((*fPIDMatrix)(index,1)); | |
711 | Int_t isPion = TMath::Nint((*fPIDMatrix)(index,2)); | |
712 | Int_t isKaon = TMath::Nint((*fPIDMatrix)(index,3)); | |
713 | Int_t isProton = TMath::Nint((*fPIDMatrix)(index,4)); | |
714 | Float_t dca[2]; | |
715 | track->GetImpactParameters(dca[0],dca[1]); | |
716 | // | |
717 | if ( (isMuon==0 && isElectron==0) && (TMath::Sqrt(dca[0]*dca[0]+dca[1]*dca[1])>kDCAcut) ) return; | |
718 | Double_t normdEdx= track->GetTPCsignal()/(medianMIP0); // TPC signal normalized to the MIP | |
719 | // | |
720 | AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam(); | |
721 | AliExternalTrackParam * trackOut = (AliExternalTrackParam *)track->GetOuterParam(); | |
722 | AliExternalTrackParam * tpcOut = (AliExternalTrackParam *)ftrack->GetTPCOut(); | |
723 | if (!trackIn) return; | |
724 | if (!trackOut) return; | |
725 | if (!tpcOut) return; | |
726 | if (trackIn->GetZ()*trackOut->GetZ()<0) return; // remove crossing tracks | |
727 | // | |
728 | // calculate local and global angle | |
729 | Int_t side = (trackIn->GetZ()>0)? 1:-1; | |
730 | Double_t tgl=trackIn->GetTgl(); | |
731 | Double_t gangle[3]={0,0,0}; | |
732 | Double_t langle[3]={0,0,0}; | |
733 | Double_t length[3]={0,0,0}; | |
734 | Double_t pxpypz[3]={0,0,0}; | |
735 | Double_t bz=fMagF; | |
736 | trackIn->GetXYZAt(radius0,bz,pxpypz); // get the global position at the middle of the IROC | |
737 | gangle[0]=TMath::ATan2(pxpypz[1],pxpypz[0]); // global angle IROC | |
738 | trackIn->GetXYZAt(radius1,bz,pxpypz); // get the global position at the middle of the OROC - medium pads | |
739 | gangle[1]=TMath::ATan2(pxpypz[1],pxpypz[0]); // global angle OROC | |
740 | trackOut->GetXYZAt(radius2,bz,pxpypz); // get the global position at the middle of OROC - long pads | |
741 | gangle[2]=TMath::ATan2(pxpypz[1],pxpypz[0]); | |
742 | // | |
743 | trackIn->GetPxPyPzAt(radius0,bz,pxpypz); //get momentum vector | |
744 | langle[0]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[0]; //local angle between padrow and track IROC | |
745 | trackIn->GetPxPyPzAt(radius1,bz,pxpypz); | |
746 | langle[1]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[1]; | |
747 | trackOut->GetPxPyPzAt(radius2,bz,pxpypz); // OROC medium | |
748 | langle[2]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[2]; | |
749 | for (Int_t i=0; i<3; i++){ | |
750 | if (langle[i]>TMath::Pi()) langle[i]-=TMath::TwoPi(); | |
751 | if (langle[i]<-TMath::Pi()) langle[i]+=TMath::TwoPi(); | |
752 | length[i]=TMath::Sqrt(1+langle[i]*langle[i]+tgl*tgl); // the tracklet length | |
753 | } | |
754 | // | |
755 | // Select the kaons and Protons which are "isolated" in TPC dedx curve | |
756 | // | |
757 | // | |
758 | Double_t dedxP = AliExternalTrackParam::BetheBlochAleph(track->GetInnerParam()->GetP()/massP,kp1,kp2,kp3,kp4,kp5); | |
759 | Double_t dedxK = AliExternalTrackParam::BetheBlochAleph(track->GetInnerParam()->GetP()/massK,kp1,kp2,kp3,kp4,kp5); | |
760 | if (dedxP>2 || dedxK>2){ | |
761 | if (track->GetP()<1.2 && normdEdx>1.8&&counterMIP0>10){ // not enough from TOF and V0s triggered by high dedx | |
762 | // signing the Proton and kaon - using the "bitmask" bit 1 and 2 is dedicated for V0s and TOF selected | |
763 | if ( TMath::Abs(normdEdx/dedxP-1)<0.3) isProton+=4; | |
764 | if ( TMath::Abs(normdEdx/dedxK-1)<0.3) isKaon+=4; | |
765 | if (normdEdx/dedxK>1.3) isProton+=8; | |
766 | if (normdEdx/dedxP<0.7) isKaon+=8; | |
767 | } | |
768 | } | |
769 | // | |
770 | // | |
771 | // | |
772 | Double_t mass = 0; | |
773 | Bool_t isHighPt = ((TMath::Power(1/track->Pt(),4)*gRandom->Rndm())<0.005); // rnadomly selected HPT tracks | |
774 | // there are selected for the QA of the obtained calibration | |
775 | Bool_t isMIP = TMath::Abs(track->GetInnerParam()->P()-0.4)<0.005&&(counter<kMax); // | |
776 | // REMINDER - it is not exactly MIP - we select the regtion where the Kaon and Electrons are well separated | |
777 | ||
778 | if (isElectron>0) mass = massE; | |
779 | if (isProton>0) mass = massP; | |
780 | if (isKaon>0) mass = massK; | |
781 | if (isMuon>0) mass = massMuon; | |
782 | if (isPion>0) mass = massPi; | |
783 | if (isHighPt) mass = massPi; //assign mass of pions | |
784 | if (isMIP&&track->GetTPCsignal()/medianMIP0<1.5) mass = massPi; //assign mass of pions | |
785 | if (mass==0) return; | |
786 | // | |
787 | // calculate expected dEdx | |
788 | Double_t dedxDef= 0; | |
789 | Double_t dedxDefPion= 0,dedxDefProton=0, dedxDefKaon=0; | |
790 | Double_t pin=trackIn->GetP(); | |
791 | Double_t pout=trackOut->GetP(); | |
792 | Double_t p=(pin+pout)*0.5; // momenta as the mean between tpc momenta at inner and outer wall of the TPC | |
793 | if (mass>0) dedxDef = AliExternalTrackParam::BetheBlochAleph(p/mass,kp1,kp2,kp3,kp4,kp5); | |
794 | dedxDefPion = AliExternalTrackParam::BetheBlochAleph(p/massPi,kp1,kp2,kp3,kp4,kp5); | |
795 | dedxDefProton = AliExternalTrackParam::BetheBlochAleph(p/massP,kp1,kp2,kp3,kp4,kp5); | |
796 | dedxDefKaon = AliExternalTrackParam::BetheBlochAleph(p/massK,kp1,kp2,kp3,kp4,kp5); | |
797 | // | |
798 | // dEdx Truncated mean vectros with differnt tuncation | |
799 | // 11 truncations array - 0-10 - 0~50% 11=100% | |
800 | // 3 Regions - IROC,OROC0, OROC1 | |
801 | // 2 Q - total charge and max charge | |
802 | // Log - Logarithmic mean used | |
803 | // Up/Dwon - Upper half or lower half of truncation used | |
804 | // RMS - rms of the distribction (otherwise truncated mean) | |
805 | // M2 suffix - second moment ( truncated) | |
806 | TVectorF truncUp(11); | |
807 | TVectorF truncDown(11); | |
808 | TVectorF vecAllMax(11); | |
809 | TVectorF vecIROCMax(11); | |
810 | TVectorF vecOROCMax(11); | |
811 | TVectorF vecOROC0Max(11); | |
812 | TVectorF vecOROC1Max(11); | |
813 | // | |
814 | TVectorF vecAllTot(11); | |
815 | TVectorF vecIROCTot(11); | |
816 | TVectorF vecOROCTot(11); | |
817 | TVectorF vecOROC0Tot(11); | |
818 | TVectorF vecOROC1Tot(11); | |
819 | // | |
820 | TVectorF vecAllTotLog(11); | |
821 | TVectorF vecIROCTotLog(11); | |
822 | TVectorF vecOROCTotLog(11); | |
823 | TVectorF vecOROC0TotLog(11); | |
824 | TVectorF vecOROC1TotLog(11); | |
825 | // | |
826 | TVectorF vecAllTotUp(11); | |
827 | TVectorF vecIROCTotUp(11); | |
828 | TVectorF vecOROCTotUp(11); | |
829 | TVectorF vecOROC0TotUp(11); | |
830 | TVectorF vecOROC1TotUp(11); | |
831 | // | |
832 | TVectorF vecAllTotDown(11); | |
833 | TVectorF vecIROCTotDown(11); | |
834 | TVectorF vecOROCTotDown(11); | |
835 | TVectorF vecOROC0TotDown(11); | |
836 | TVectorF vecOROC1TotDown(11); | |
837 | ||
838 | TVectorF vecAllTotRMS(11); | |
839 | TVectorF vecIROCTotRMS(11); | |
840 | TVectorF vecOROCTotRMS(11); | |
841 | TVectorF vecOROC0TotRMS(11); | |
842 | TVectorF vecOROC1TotRMS(11); | |
843 | // | |
844 | TVectorF vecAllTotM2(11); | |
845 | TVectorF vecIROCTotM2(11); | |
846 | TVectorF vecOROCTotM2(11); | |
847 | TVectorF vecOROC0TotM2(11); | |
848 | TVectorF vecOROC1TotM2(11); | |
849 | // | |
850 | TVectorF vecAllTotMS(11); | |
851 | TVectorF vecIROCTotMS(11); | |
852 | TVectorF vecOROCTotMS(11); | |
853 | TVectorF vecOROC0TotMS(11); | |
854 | TVectorF vecOROC1TotMS(11); | |
855 | // | |
856 | // Differnt number of clusters definitions - in separate regions of the TPC | |
857 | // 20 - ratio - found/findabel | |
858 | // 21 - number of clusters used for given dEdx calculation | |
859 | // | |
860 | // suffix - 3 or 4 - number of padrows before and after given row to define findable row | |
861 | // | |
862 | Double_t ncl20All = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3); | |
863 | Double_t ncl20IROC = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3); | |
864 | Double_t ncl20OROC = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3); | |
865 | Double_t ncl20OROC0= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3); | |
866 | Double_t ncl20OROC1= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3); | |
867 | // | |
868 | Double_t ncl20All4 = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3,4); | |
869 | Double_t ncl20IROC4 = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3,4); | |
870 | Double_t ncl20OROC4 = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3,4); | |
871 | Double_t ncl20OROC04= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3,4); | |
872 | Double_t ncl20OROC14= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3,4); | |
873 | // | |
874 | Double_t ncl20All3 = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3,3); | |
875 | Double_t ncl20IROC3 = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3,3); | |
876 | Double_t ncl20OROC3 = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3,3); | |
877 | Double_t ncl20OROC03= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3,3); | |
878 | Double_t ncl20OROC13= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3,3); | |
879 | // | |
880 | Double_t ncl21All = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,2); | |
881 | Double_t ncl21IROC = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,2); | |
882 | Double_t ncl21OROC = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,2); | |
883 | Double_t ncl21OROC0= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,2); | |
884 | Double_t ncl21OROC1= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,2); | |
885 | // calculate truncated dEdx - mean rms M2 ... | |
886 | Int_t ifrac=0; | |
887 | for (Int_t ifracDown=0; ifracDown<1; ifracDown++){ | |
888 | for (Int_t ifracUp=0; ifracUp<11; ifracUp++){ | |
889 | Double_t fracDown = 0.0+Double_t(ifracDown)*0.05; | |
890 | Double_t fracUp = 0.5+Double_t(ifracUp)*0.05; | |
891 | vecAllMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,0,159,0); | |
892 | vecIROCMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,0,63,0); | |
893 | vecOROCMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,64,159,0); | |
894 | vecOROC0Max[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,64,128,0); | |
895 | vecOROC1Max[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,129,159,0); | |
896 | // | |
897 | vecAllTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,0); | |
898 | vecIROCTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,0); | |
899 | vecOROCTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,0); | |
900 | vecOROC0Tot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,0); | |
901 | vecOROC1Tot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,0); | |
902 | // | |
903 | vecAllTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,0,2,1); | |
904 | vecIROCTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,0,2,1); | |
905 | vecOROCTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,0,2,1); | |
906 | vecOROC0TotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,0,2,1); | |
907 | vecOROC1TotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,0,2,1); | |
908 | // | |
909 | vecAllTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,4,2,1); | |
910 | vecIROCTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,4,2,1); | |
911 | vecOROCTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,4,2,1); | |
912 | vecOROC0TotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,4,2,1); | |
913 | vecOROC1TotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,4,2,1); | |
914 | // | |
915 | vecAllTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,5,2,1); | |
916 | vecIROCTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,5,2,1); | |
917 | vecOROCTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,5,2,1); | |
918 | vecOROC0TotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,5,2,1); | |
919 | vecOROC1TotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,5,2,1); | |
920 | // | |
921 | vecAllTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,1,2,0); | |
922 | vecIROCTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,1,2,0); | |
923 | vecOROCTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,1,2,0); | |
924 | vecOROC0TotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,1,2,0); | |
925 | vecOROC1TotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,1,2,0); | |
926 | // | |
927 | vecAllTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,6,2,1); | |
928 | vecIROCTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,6,2,1); | |
929 | vecOROCTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,6,2,1); | |
930 | vecOROC0TotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,6,2,1); | |
931 | vecOROC1TotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,6,2,1); | |
932 | // | |
933 | vecAllTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,8,2,1); | |
934 | vecIROCTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,8,2,1); | |
935 | vecOROCTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,8,2,1); | |
936 | vecOROC0TotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,8,2,1); | |
937 | vecOROC1TotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,8,2,1); | |
938 | truncUp[ifrac]=fracUp; | |
939 | truncDown[ifrac]=fracDown; | |
940 | ifrac++; | |
941 | } | |
942 | } | |
943 | // | |
944 | // Fill histograms | |
945 | // | |
946 | if ((isKaon||isProton||isPion||isElectron||isMIP||isMuon&&(!isHighPt)) && dedxDef>0) { | |
947 | // | |
948 | Int_t ncont = vertex->GetNContributors(); | |
949 | for (Int_t ipad=0; ipad<3; ipad++){ | |
950 | // histogram with enahanced phi granularity - to make gain phi maps | |
951 | Double_t xxx[4]={0,gangle[ipad],1./dedxDef,ipad*2+((side>0)?0:1)}; | |
952 | Double_t nclR=0; | |
953 | if (ipad==0) {xxx[0]=vecIROCTot[4]/medianMIP0; nclR=ncl21IROC/63.;} | |
954 | if (ipad==1) {xxx[0]=vecOROC0Tot[4]/medianMIP0;nclR=ncl21OROC0/63.;} | |
955 | if (ipad==2) {xxx[0]=vecOROC1Tot[4]/medianMIP0;nclR=ncl21OROC1/32.;} | |
956 | xxx[0]/=dedxDef; | |
957 | if (xxx[0]>0) xxx[0]=1./xxx[0]; | |
958 | if (TMath::Abs(langle[ipad])<0.25&&nclR>0.4) fHistdEdxMap->Fill(xxx); | |
959 | } | |
960 | for (Int_t ipad=0; ipad<3; ipad++){ | |
961 | // | |
962 | // this are histogram to define overall main gain correction | |
963 | // Maybe dead end - we can not put all info which can be used into the THnSparse | |
964 | // It is keeped there for educational point of view | |
965 | // | |
966 | Double_t xxx[6]={0,1./dedxDef, TMath::Abs(langle[ipad]), TMath::Abs(tgl), ncont, ipad }; | |
967 | if (ipad==0) {xxx[0]=vecIROCTot[4]/medianMIP0;} | |
968 | if (ipad==1) {xxx[0]=vecOROC0Tot[4]/medianMIP0;} | |
969 | if (ipad==2) {xxx[0]=vecOROC1Tot[4]/medianMIP0;} | |
970 | xxx[0]/=dedxDef; | |
971 | if (xxx[0]>0) xxx[0]=1./xxx[0]; | |
972 | if (xxx[0]>0) fHistdEdxTot->Fill(xxx); | |
973 | if (ipad==0) {xxx[0]=vecIROCMax[4]/medianMIP0;} | |
974 | if (ipad==1) {xxx[0]=vecOROC0Max[4]/medianMIP0;} | |
975 | if (ipad==2) {xxx[0]=vecOROC1Max[4]/medianMIP0;} | |
976 | xxx[0]=dedxDef; | |
977 | if (xxx[0]>0) xxx[0]=1./xxx[0]; | |
978 | fHistdEdxMax->Fill(xxx); | |
979 | } | |
980 | } | |
981 | // | |
982 | // Downscale selected tracks before filling the tree | |
983 | // | |
984 | Bool_t isSelected = kFALSE; | |
985 | // | |
986 | if (isKaon||isProton||isPion||isElectron||isMIP||isMuon) isSelected=kTRUE; | |
987 | isHighPt = kFALSE; | |
988 | if (!isSelected) isHighPt = ((TMath::Power(1/track->Pt(),4)*gRandom->Rndm())<0.005); | |
989 | if (counter>kMax && ((1/track->Pt()*gRandom->Rndm())>kMax/counter)) return; | |
990 | isSelected|=isHighPt; | |
991 | // | |
992 | // | |
993 | // | |
994 | // Equalize statistic in BB bins - special care about pions | |
995 | Int_t entriesBB = (Int_t)hisBB->GetEntries(); | |
996 | if ((isElectron==0 &&isMuon==0 && p<2.) && entriesBB>20 &&dedxDef>0.01){ | |
997 | Int_t bin = hisBB->GetXaxis()->FindBin(1./dedxDef); | |
998 | Double_t cont = hisBB->GetBinContent(bin); | |
999 | Double_t mean =(entriesBB)/20.; | |
1000 | if ((isPion>0) && gRandom->Rndm()*cont > 0.1*mean) return; | |
1001 | if ((isPion==0) && gRandom->Rndm()*cont > 0.25*mean) return; | |
1002 | } | |
1003 | if (!isSelected) return; | |
1004 | if (dedxDef>0.01) hisBB->Fill(1./dedxDef); | |
1005 | // | |
1006 | if (isHighPt) counterHPT++; | |
1007 | counter++; | |
1008 | // | |
1009 | TTreeSRedirector * pcstream = GetDebugStreamer(); | |
1010 | Double_t ptrel0 = AliTPCcalibDB::GetPTRelative(fTime,fRun,0); | |
1011 | Double_t ptrel1 = AliTPCcalibDB::GetPTRelative(fTime,fRun,1); | |
1012 | Int_t sectorIn = Int_t(18+9*(trackIn->GetAlpha()/TMath::Pi()))%18; | |
1013 | Int_t sectorOut = Int_t(18+9*(trackOut->GetAlpha()/TMath::Pi()))%18; | |
1014 | // | |
1015 | if (pcstream){ | |
1016 | (*pcstream)<<"dump"<< | |
1017 | "vertex.="<<vertex<< | |
1018 | "bz="<<fMagF<< | |
1019 | "ptrel0="<<ptrel0<< | |
1020 | "ptrel1="<<ptrel1<< | |
1021 | "sectorIn="<<sectorIn<< | |
1022 | "sectorOut="<<sectorOut<< | |
1023 | "side="<<side<< | |
1024 | // pid type | |
1025 | "isMuon="<<isMuon<< | |
1026 | "isProton="<<isProton<< | |
1027 | "isKaon="<<isKaon<< | |
1028 | "isPion="<<isPion<< | |
1029 | "isElectron="<<isElectron<< | |
1030 | "isMIP="<<isMIP<< | |
1031 | "isHighPt="<<isHighPt<< | |
1032 | "mass="<<mass<< | |
1033 | "dedxDef="<<dedxDef<< | |
1034 | "dedxDefPion="<<dedxDefPion<< | |
1035 | "dedxDefKaon="<<dedxDefKaon<< | |
1036 | "dedxDefProton="<<dedxDefProton<< | |
1037 | // | |
1038 | "nclITS="<<nclITS<< | |
1039 | "ncl="<<ncl<< | |
1040 | "ncl21="<<ncl21<< | |
1041 | "ncl20="<<ncl20<< | |
1042 | // | |
1043 | "ncl20All="<<ncl20All<< | |
1044 | "ncl20OROC="<<ncl20OROC<< | |
1045 | "ncl20IROC="<<ncl20IROC<< | |
1046 | "ncl20OROC0="<<ncl20OROC0<< | |
1047 | "ncl20OROC1="<<ncl20OROC1<< | |
1048 | // | |
1049 | "ncl20All4="<<ncl20All4<< | |
1050 | "ncl20OROC4="<<ncl20OROC4<< | |
1051 | "ncl20IROC4="<<ncl20IROC4<< | |
1052 | "ncl20OROC04="<<ncl20OROC04<< | |
1053 | "ncl20OROC14="<<ncl20OROC14<< | |
1054 | // | |
1055 | "ncl20All3="<<ncl20All3<< | |
1056 | "ncl20OROC3="<<ncl20OROC3<< | |
1057 | "ncl20IROC3="<<ncl20IROC3<< | |
1058 | "ncl20OROC03="<<ncl20OROC03<< | |
1059 | "ncl20OROC13="<<ncl20OROC13<< | |
1060 | // | |
1061 | "ncl21All="<<ncl21All<< | |
1062 | "ncl21OROC="<<ncl21OROC<< | |
1063 | "ncl21IROC="<<ncl21IROC<< | |
1064 | "ncl21OROC0="<<ncl21OROC0<< | |
1065 | "ncl21OROC1="<<ncl21OROC1<< | |
1066 | //track properties | |
1067 | "langle0="<<langle[0]<< | |
1068 | "langle1="<<langle[1]<< | |
1069 | "langle2="<<langle[2]<< | |
1070 | "gangle0="<<gangle[0]<< //global angle phi IROC | |
1071 | "gangle1="<<gangle[1]<< // OROC medium | |
1072 | "gangle2="<<gangle[2]<< // OROC long | |
1073 | "L0="<<length[0]<< | |
1074 | "L1="<<length[1]<< | |
1075 | "L2="<<length[2]<< | |
1076 | "p="<<p<< | |
1077 | "pin="<<pin<< | |
1078 | "pout="<<pout<< | |
1079 | "tgl="<<tgl<< | |
1080 | "track.="<<track<< | |
1081 | "trackIn.="<<trackIn<< | |
1082 | "trackOut.="<<trackOut<< | |
1083 | "tpcOut.="<<tpcOut<< | |
1084 | "medianMIP0="<<medianMIP0<< // median MIP position as estimated from the array of (not only) "MIPS" | |
1085 | //dedx | |
1086 | "truncUp.="<<&truncUp<< | |
1087 | "truncDown.="<<&truncDown<< | |
1088 | "vecAllMax.="<<&vecAllMax<< | |
1089 | "vecIROCMax.="<<&vecIROCMax<< | |
1090 | "vecOROCMax.="<<&vecOROCMax<< | |
1091 | "vecOROC0Max.="<<&vecOROC0Max<< | |
1092 | "vecOROC1Max.="<<&vecOROC1Max<< | |
1093 | // | |
1094 | "vecAllTot.="<<&vecAllTot<< | |
1095 | "vecIROCTot.="<<&vecIROCTot<< | |
1096 | "vecOROCTot.="<<&vecOROCTot<< | |
1097 | "vecOROC0Tot.="<<&vecOROC0Tot<< | |
1098 | "vecOROC1Tot.="<<&vecOROC1Tot<< | |
1099 | // | |
1100 | "vecAllTotLog.="<<&vecAllTotLog<< | |
1101 | "vecIROCTotLog.="<<&vecIROCTotLog<< | |
1102 | "vecOROCTotLog.="<<&vecOROCTotLog<< | |
1103 | "vecOROC0TotLog.="<<&vecOROC0TotLog<< | |
1104 | "vecOROC1TotLog.="<<&vecOROC1TotLog<< | |
1105 | // | |
1106 | "vecAllTotUp.="<<&vecAllTotUp<< | |
1107 | "vecIROCTotUp.="<<&vecIROCTotUp<< | |
1108 | "vecOROCTotUp.="<<&vecOROCTotUp<< | |
1109 | "vecOROC0TotUp.="<<&vecOROC0TotUp<< | |
1110 | "vecOROC1TotUp.="<<&vecOROC1TotUp<< | |
1111 | // | |
1112 | "vecAllTotDown.="<<&vecAllTotDown<< | |
1113 | "vecIROCTotDown.="<<&vecIROCTotDown<< | |
1114 | "vecOROCTotDown.="<<&vecOROCTotDown<< | |
1115 | "vecOROC0TotDown.="<<&vecOROC0TotDown<< | |
1116 | "vecOROC1TotDown.="<<&vecOROC1TotDown<< | |
1117 | // | |
1118 | "vecAllTotRMS.="<<&vecAllTotRMS<< | |
1119 | "vecIROCTotRMS.="<<&vecIROCTotRMS<< | |
1120 | "vecOROCTotRMS.="<<&vecOROCTotRMS<< | |
1121 | "vecOROC0TotRMS.="<<&vecOROC0TotRMS<< | |
1122 | "vecOROC1TotRMS.="<<&vecOROC1TotRMS<< | |
1123 | // | |
1124 | "vecAllTotM2.="<<&vecAllTotM2<< | |
1125 | "vecIROCTotM2.="<<&vecIROCTotM2<< | |
1126 | "vecOROCTotM2.="<<&vecOROCTotM2<< | |
1127 | "vecOROC0TotM2.="<<&vecOROC0TotM2<< | |
1128 | "vecOROC1TotM2.="<<&vecOROC1TotM2<< | |
1129 | // | |
1130 | "vecAllTotMS.="<<&vecAllTotMS<< | |
1131 | "vecIROCTotMS.="<<&vecIROCTotMS<< | |
1132 | "vecOROCTotMS.="<<&vecOROCTotMS<< | |
1133 | "vecOROC0TotMS.="<<&vecOROC0TotMS<< | |
1134 | "vecOROC1TotMS.="<<&vecOROC1TotMS<< | |
1135 | "\n"; | |
1136 | } | |
1137 | } | |
1138 | ||
1139 | ||
1140 | ||
1141 | ||
1142 | void AliTPCcalibGainMult::ProcessV0s(AliESDEvent * event){ | |
1143 | // | |
1144 | // Select the K0s and gamma - and sign daughter products | |
1145 | // | |
1146 | TTreeSRedirector * pcstream = GetDebugStreamer(); | |
1147 | AliKFParticle::SetField(event->GetMagneticField()); | |
1148 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend")); | |
1149 | if (!esdFriend) { | |
1150 | //Printf("ERROR: esdFriend not available"); | |
1151 | return; | |
1152 | } | |
1153 | if (esdFriend->TestSkipBit()) return; | |
1154 | // | |
1155 | // | |
1156 | TDatabasePDG *pdg = TDatabasePDG::Instance(); | |
1157 | const Double_t kChi2Cut=5; | |
1158 | const Double_t kMinR=2; | |
1159 | const Int_t kMinNcl=110; | |
1160 | const Double_t kMinREl=5; | |
1161 | const Double_t kMaxREl=70; | |
1162 | // | |
1163 | Int_t nv0 = event->GetNumberOfV0s(); | |
1164 | AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex(); | |
1165 | AliKFVertex kfvertex=*vertex; | |
1166 | // | |
1167 | for (Int_t iv0=0;iv0<nv0;iv0++){ | |
1168 | AliESDv0 *v0 = event->GetV0(iv0); | |
1169 | if (!v0) continue; | |
1170 | if (v0->GetOnFlyStatus()<0.5) continue; | |
1171 | if (v0->GetPindex()<0) continue; | |
1172 | if (v0->GetNindex()<0) continue; | |
1173 | if (TMath::Max(v0->GetPindex(), v0->GetNindex())>event->GetNumberOfTracks()) continue; | |
1174 | // | |
1175 | // | |
1176 | AliExternalTrackParam pp=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamP())):(*(v0->GetParamN())); | |
1177 | AliExternalTrackParam pn=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamN())):(*(v0->GetParamP())); | |
1178 | AliKFParticle kfp1( pp, 211 ); | |
1179 | AliKFParticle kfp2( pn, -211 ); | |
1180 | // | |
1181 | AliKFParticle *v0KFK0 = new AliKFParticle(kfp1,kfp2); | |
1182 | AliKFParticle *v0KFK0CV = new AliKFParticle(*v0KFK0); | |
1183 | v0KFK0CV->SetProductionVertex(kfvertex); | |
1184 | v0KFK0CV->TransportToProductionVertex(); | |
1185 | Double_t chi2K0 = v0KFK0CV->GetChi2(); | |
1186 | if (chi2K0>kChi2Cut) continue; | |
1187 | if (v0->GetRr()<kMinR) continue; | |
1188 | Bool_t isOKC=TMath::Max(v0->GetCausalityP()[0],v0->GetCausalityP()[1])<0.7&&TMath::Min(v0->GetCausalityP()[2],v0->GetCausalityP()[3])>0.2; | |
1189 | // | |
1190 | Double_t effMass22=v0->GetEffMass(2,2); | |
1191 | Double_t effMass42=v0->GetEffMass(4,2); | |
1192 | Double_t effMass24=v0->GetEffMass(2,4); | |
1193 | Double_t effMass00=v0->GetEffMass(0,0); | |
1194 | AliKFParticle *v0KFK0CVM = new AliKFParticle(*v0KFK0CV); | |
1195 | v0KFK0CVM->SetMassConstraint(pdg->GetParticle("K_S0")->Mass()); | |
1196 | Bool_t isV0= kFALSE; | |
1197 | // | |
1198 | Double_t d22 = TMath::Abs(effMass22-pdg->GetParticle("K_S0")->Mass()); | |
1199 | Double_t d42 = TMath::Abs(effMass42-pdg->GetParticle("Lambda0")->Mass()); | |
1200 | Double_t d24 = TMath::Abs(effMass24-pdg->GetParticle("Lambda0")->Mass()); | |
1201 | Double_t d00 = TMath::Abs(effMass00); | |
1202 | // | |
1203 | Bool_t isKaon = d22<0.01 && d22< 0.3 * TMath::Min(TMath::Min(d42,d24),d00); | |
1204 | Bool_t isLambda = d42<0.01 && d42< 0.3 * TMath::Min(TMath::Min(d22,d24),d00); | |
1205 | Bool_t isAntiLambda= d24<0.01 && d24< 0.3 * TMath::Min(TMath::Min(d22,d42),d00); | |
1206 | Bool_t isGamma = d00<0.02 && d00< 0.3 * TMath::Min(TMath::Min(d42,d24),d22); | |
1207 | // | |
1208 | if (isGamma && (isKaon||isLambda||isAntiLambda)) continue; | |
1209 | if (isLambda && (isKaon||isGamma||isAntiLambda)) continue; | |
1210 | if (isKaon && (isLambda||isGamma||isAntiLambda)) continue; | |
1211 | Double_t sign= v0->GetParamP()->GetSign()* v0->GetParamN()->GetSign(); | |
1212 | if (sign>0) continue; | |
1213 | isV0=isKaon||isLambda||isAntiLambda||isGamma; | |
1214 | if (!(isV0)) continue; | |
1215 | if (isGamma&&v0->GetRr()<kMinREl) continue; | |
1216 | if (isGamma&&v0->GetRr()>kMaxREl) continue; | |
1217 | if (!isOKC) continue; | |
1218 | // | |
1219 | Int_t pindex = (v0->GetParamP()->GetSign()>0) ? v0->GetPindex() : v0->GetNindex(); | |
1220 | Int_t nindex = (v0->GetParamP()->GetSign()>0) ? v0->GetNindex() : v0->GetPindex(); | |
1221 | AliESDtrack * trackP = event->GetTrack(pindex); | |
1222 | AliESDtrack * trackN = event->GetTrack(nindex); | |
1223 | if (!trackN) continue; | |
1224 | if (!trackP) continue; | |
1225 | Int_t nclP= (Int_t)trackP->GetTPCClusterInfo(2,1); | |
1226 | Int_t nclN= (Int_t)trackN->GetTPCClusterInfo(2,1); | |
1227 | if (TMath::Min(nclP,nclN)<kMinNcl) continue; | |
1228 | Double_t eta = TMath::Max(TMath::Abs(trackP->Eta()), TMath::Abs(trackN->Eta())); | |
1229 | if (TMath::Abs(eta)>1) continue; | |
1230 | // | |
1231 | // | |
1232 | AliESDfriendTrack *friendTrackP = esdFriend->GetTrack(pindex); | |
1233 | AliESDfriendTrack *friendTrackN = esdFriend->GetTrack(nindex); | |
1234 | if (!friendTrackP) continue; | |
1235 | if (!friendTrackN) continue; | |
1236 | TObject *calibObject; | |
1237 | AliTPCseed *seedP = 0; | |
1238 | AliTPCseed *seedN = 0; | |
1239 | for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) { | |
1240 | if ((seedP=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
1241 | } | |
1242 | for (Int_t l=0;(calibObject=friendTrackN->GetCalibObject(l));++l) { | |
1243 | if ((seedN=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
1244 | } | |
1245 | if (isGamma){ | |
1246 | if ( TMath::Abs((trackP->GetTPCsignal()/(trackN->GetTPCsignal()+0.0001)-1)>0.3)) continue; | |
1247 | } | |
1248 | if (isGamma) (*fPIDMatrix)(pindex, 0)+=2; | |
1249 | if (isGamma) (*fPIDMatrix)(nindex, 0)+=2; | |
1250 | // | |
1251 | if (isKaon) (*fPIDMatrix)(pindex, 2)+=2; | |
1252 | if (isKaon) (*fPIDMatrix)(nindex, 2)+=2; | |
1253 | // | |
1254 | // | |
1255 | if (pcstream){ | |
1256 | (*pcstream)<<"v0s"<< | |
1257 | "isGamma="<<isGamma<< | |
1258 | "isKaon="<<isKaon<< | |
1259 | "isLambda="<<isLambda<< | |
1260 | "isAntiLambda="<<isAntiLambda<< | |
1261 | "chi2="<<chi2K0<< | |
1262 | "trackP.="<<trackP<< | |
1263 | "trackN.="<<trackN<< | |
1264 | "v0.="<<v0<< | |
1265 | "\n"; | |
1266 | } | |
1267 | } | |
1268 | } | |
1269 | ||
1270 | ||
1271 | ||
1272 | ||
1273 | void AliTPCcalibGainMult::ProcessCosmic(const AliESDEvent * event) { | |
1274 | // | |
1275 | // Find cosmic pairs trigger by random trigger | |
1276 | // | |
1277 | // | |
1278 | AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; | |
1279 | AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); | |
1280 | ||
1281 | AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD(); | |
1282 | AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC(); | |
1283 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend")); | |
1284 | const Double_t kMinPt=4; | |
1285 | const Double_t kMinPtMax=0.8; | |
1286 | const Double_t kMinNcl=159*0.5; | |
1287 | const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1}; | |
1288 | Int_t ntracks=event->GetNumberOfTracks(); | |
1289 | const Double_t kMaxImpact=80; | |
1290 | // Float_t dcaTPC[2]={0,0}; | |
1291 | // Float_t covTPC[3]={0,0,0}; | |
1292 | ||
1293 | UInt_t specie = event->GetEventSpecie(); // skip laser events | |
1294 | if (specie==AliRecoParam::kCalib) return; | |
1295 | ||
1296 | ||
1297 | for (Int_t itrack0=0;itrack0<ntracks;itrack0++) { | |
1298 | AliESDtrack *track0 = event->GetTrack(itrack0); | |
1299 | if (!track0) continue; | |
1300 | if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue; | |
1301 | ||
1302 | if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()<kMinPt) continue; | |
1303 | if (track0->GetTPCncls()<kMinNcl) continue; | |
1304 | if (TMath::Abs(track0->GetY())<2*kMaxDelta[0]) continue; | |
1305 | if (TMath::Abs(track0->GetY())>kMaxImpact) continue; | |
1306 | if (track0->GetKinkIndex(0)>0) continue; | |
1307 | const Double_t * par0=track0->GetParameter(); //track param at rhe DCA | |
1308 | //rm primaries | |
1309 | // | |
1310 | for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) { | |
1311 | AliESDtrack *track1 = event->GetTrack(itrack1); | |
1312 | if (!track1) continue; | |
1313 | if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue; | |
1314 | if (track1->GetKinkIndex(0)>0) continue; | |
1315 | if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()<kMinPt) continue; | |
1316 | if (track1->GetTPCncls()<kMinNcl) continue; | |
1317 | if (TMath::Abs(AliTracker::GetBz())>1&&TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue; | |
1318 | if (TMath::Abs(track1->GetY())<2*kMaxDelta[0]) continue; | |
1319 | if (TMath::Abs(track1->GetY())>kMaxImpact) continue; | |
1320 | // | |
1321 | const Double_t* par1=track1->GetParameter(); //track param at rhe DCA | |
1322 | // | |
1323 | Bool_t isPair=kTRUE; | |
1324 | for (Int_t ipar=0; ipar<5; ipar++){ | |
1325 | if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off | |
1326 | if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE; | |
1327 | } | |
1328 | if (!isPair) continue; | |
1329 | if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE; | |
1330 | //delta with correct sign | |
1331 | if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign | |
1332 | if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign | |
1333 | if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign | |
1334 | if (!isPair) continue; | |
1335 | TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); | |
1336 | Int_t eventNumber = event->GetEventNumberInFile(); | |
1337 | Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(itrack0)!=0):0; | |
1338 | Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4); | |
1339 | printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS); | |
1340 | // | |
1341 | // | |
1342 | TTreeSRedirector * pcstream = GetDebugStreamer(); | |
1343 | Int_t ntracksSPD = vertexSPD->GetNContributors(); | |
1344 | Int_t ntracksTPC = vertexTPC->GetNContributors(); | |
1345 | // | |
1346 | if (pcstream){ | |
1347 | (*pcstream)<<"cosmicPairsAll"<< | |
1348 | "run="<<fRun<< // run number | |
1349 | "event="<<fEvent<< // event number | |
1350 | "time="<<fTime<< // time stamp of event | |
1351 | "trigger="<<fTrigger<< // trigger | |
1352 | "triggerClass="<<&fTriggerClass<< // trigger | |
1353 | "bz="<<fMagF<< // magnetic field | |
1354 | // | |
1355 | "nSPD="<<ntracksSPD<< | |
1356 | "nTPC="<<ntracksTPC<< | |
1357 | "vSPD.="<<vertexSPD<< //primary vertex -SPD | |
1358 | "vTPC.="<<vertexTPC<< //primary vertex -TPC | |
1359 | "t0.="<<track0<< //track0 | |
1360 | "t1.="<<track1<< //track1 | |
1361 | "\n"; | |
1362 | } | |
1363 | // | |
1364 | AliESDfriendTrack *friendTrack0 = esdFriend->GetTrack(itrack0); | |
1365 | if (!friendTrack0) continue; | |
1366 | AliESDfriendTrack *friendTrack1 = esdFriend->GetTrack(itrack1); | |
1367 | if (!friendTrack1) continue; | |
1368 | TObject *calibObject; | |
1369 | AliTPCseed *seed0 = 0; | |
1370 | AliTPCseed *seed1 = 0; | |
1371 | // | |
1372 | for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) { | |
1373 | if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
1374 | } | |
1375 | for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) { | |
1376 | if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
1377 | } | |
1378 | // | |
1379 | if (pcstream){ | |
1380 | (*pcstream)<<"cosmicPairs"<< | |
1381 | "run="<<fRun<< // run number | |
1382 | "event="<<fEvent<< // event number | |
1383 | "time="<<fTime<< // time stamp of event | |
1384 | "trigger="<<fTrigger<< // trigger | |
1385 | "triggerClass="<<&fTriggerClass<< // trigger | |
1386 | "bz="<<fMagF<< // magnetic field | |
1387 | // | |
1388 | "nSPD="<<ntracksSPD<< | |
1389 | "nTPC="<<ntracksTPC<< | |
1390 | "vSPD.="<<vertexSPD<< //primary vertex -SPD | |
1391 | "vTPC.="<<vertexTPC<< //primary vertex -TPC | |
1392 | "t0.="<<track0<< //track0 | |
1393 | "t1.="<<track1<< //track1 | |
1394 | "ft0.="<<friendTrack0<< //track0 | |
1395 | "ft1.="<<friendTrack1<< //track1 | |
1396 | "s0.="<<seed0<< //track0 | |
1397 | "s1.="<<seed1<< //track1 | |
1398 | "\n"; | |
1399 | } | |
1400 | if (!seed0) continue; | |
1401 | if (!seed1) continue; | |
1402 | Int_t nclA0=0, nclC0=0; // number of clusters | |
1403 | Int_t nclA1=0, nclC1=0; // number of clusters | |
1404 | // | |
1405 | for (Int_t irow=0; irow<159; irow++){ | |
1406 | AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow); | |
1407 | AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow); | |
1408 | if (cluster0){ | |
1409 | if (cluster0->GetQ()>0 && cluster0->GetDetector()%36<18) nclA0++; | |
1410 | if (cluster0->GetQ()>0 && cluster0->GetDetector()%36>=18) nclC0++; | |
1411 | } | |
1412 | if (cluster1){ | |
1413 | if (cluster1->GetQ()>0 && cluster1->GetDetector()%36<18) nclA1++; | |
1414 | if (cluster1->GetQ()>0 && cluster1->GetDetector()%36>=18) nclC1++; | |
1415 | } | |
1416 | } | |
1417 | Int_t cosmicType=0; // types of cosmic topology | |
1418 | if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side | |
1419 | if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side | |
1420 | if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side | |
1421 | if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side | |
1422 | if (cosmicType<2) continue; // use only crossing tracks | |
1423 | // | |
1424 | Double_t deltaTimeCluster=0; | |
1425 | deltaTimeCluster=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth(); | |
1426 | if (nclA0>nclC0) deltaTimeCluster*=-1; // if A side track | |
1427 | // | |
1428 | for (Int_t irow=0; irow<159; irow++){ | |
1429 | AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow); | |
1430 | if (cluster0 &&cluster0->GetX()>10){ | |
1431 | Double_t x0[3]={cluster0->GetRow(),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster}; | |
1432 | Int_t index0[1]={cluster0->GetDetector()}; | |
1433 | transform->Transform(x0,index0,0,1); | |
1434 | cluster0->SetX(x0[0]); | |
1435 | cluster0->SetY(x0[1]); | |
1436 | cluster0->SetZ(x0[2]); | |
1437 | // | |
1438 | } | |
1439 | AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow); | |
1440 | if (cluster1&&cluster1->GetX()>10){ | |
1441 | Double_t x1[3]={cluster1->GetRow(),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster}; | |
1442 | Int_t index1[1]={cluster1->GetDetector()}; | |
1443 | transform->Transform(x1,index1,0,1); | |
1444 | cluster1->SetX(x1[0]); | |
1445 | cluster1->SetY(x1[1]); | |
1446 | cluster1->SetZ(x1[2]); | |
1447 | } | |
1448 | } | |
1449 | // | |
1450 | // | |
1451 | if (fPIDMatrix){ | |
1452 | (*fPIDMatrix)(itrack0,1)+=4; // | |
1453 | (*fPIDMatrix)(itrack1,1)+=4; // | |
1454 | } | |
1455 | } | |
1456 | } | |
1457 | } | |
1458 | ||
1459 | ||
1460 | ||
1461 | void AliTPCcalibGainMult::ProcessKinks(const AliESDEvent * event){ | |
1462 | // | |
1463 | // | |
1464 | // | |
1465 | AliKFParticle::SetField(event->GetMagneticField()); | |
1466 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend")); | |
1467 | if (!esdFriend) { | |
1468 | //Printf("ERROR: esdFriend not available"); | |
1469 | return; | |
1470 | } | |
1471 | // if (esdFriend->TestSkipBit()) return; | |
1472 | // | |
1473 | // | |
1474 | const Double_t kChi2Cut=10; | |
1475 | const Double_t kMinR=100; | |
1476 | const Double_t kMaxR=230; | |
1477 | const Int_t kMinNcl=110; | |
1478 | // | |
1479 | Int_t nkinks = event->GetNumberOfKinks(); | |
1480 | AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex(); | |
1481 | AliKFVertex kfvertex=*vertex; | |
1482 | TTreeSRedirector * pcstream = GetDebugStreamer(); | |
1483 | // | |
1484 | for (Int_t ikink=0;ikink<nkinks;ikink++){ | |
1485 | AliESDkink *kink = event->GetKink(ikink); | |
1486 | if (!kink) continue; | |
1487 | if (kink->GetIndex(0)<0) continue; | |
1488 | if (kink->GetIndex(1)<0) continue; | |
1489 | if (TMath::Max(kink->GetIndex(1), kink->GetIndex(0))>event->GetNumberOfTracks()) continue; | |
1490 | // | |
1491 | // | |
1492 | AliExternalTrackParam pd=kink->RefParamDaughter(); | |
1493 | AliExternalTrackParam pm=kink->RefParamMother(); | |
1494 | AliKFParticle kfpd( pd, 211 ); | |
1495 | AliKFParticle kfpm( pm, -13 ); | |
1496 | // | |
1497 | AliKFParticle *v0KF = new AliKFParticle(kfpm,kfpd); | |
1498 | v0KF->SetVtxGuess(kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]); | |
1499 | Double_t chi2 = v0KF->GetChi2(); | |
1500 | AliESDtrack * trackM = event->GetTrack(kink->GetIndex(0)); | |
1501 | AliESDtrack * trackD = event->GetTrack(kink->GetIndex(1)); | |
1502 | if (!trackM) continue; | |
1503 | if (!trackD) continue; | |
1504 | Int_t nclM= (Int_t)trackM->GetTPCClusterInfo(2,1); | |
1505 | Int_t nclD= (Int_t)trackD->GetTPCClusterInfo(2,1); | |
1506 | Double_t eta = TMath::Max(TMath::Abs(trackM->Eta()), TMath::Abs(trackD->Eta())); | |
1507 | Double_t kx= v0KF->GetX(); | |
1508 | Double_t ky= v0KF->GetY(); | |
1509 | Double_t kz= v0KF->GetZ(); | |
1510 | Double_t ex= v0KF->GetErrX(); | |
1511 | Double_t ey= v0KF->GetErrY(); | |
1512 | Double_t ez= v0KF->GetErrZ(); | |
1513 | // | |
1514 | Double_t radius=TMath::Sqrt(kx*kx+ky*ky); | |
1515 | Double_t alpha=TMath::ATan2(ky,kx); | |
1516 | if (!pd.Rotate(alpha)) continue; | |
1517 | if (!pm.Rotate(alpha)) continue; | |
1518 | if (!pd.PropagateTo(radius,event->GetMagneticField())) continue; | |
1519 | if (!pm.PropagateTo(radius,event->GetMagneticField())) continue; | |
1520 | Double_t pos[2]={0,kz}; | |
1521 | Double_t cov[3]={ex*ex+ey*ey,0,ez*ez}; | |
1522 | pd.Update(pos,cov); | |
1523 | pm.Update(pos,cov); | |
1524 | // | |
1525 | if (pcstream){ | |
1526 | (*pcstream)<<"kinks"<< | |
1527 | "eta="<<eta<< | |
1528 | "nclM="<<nclM<< | |
1529 | "nclD="<<nclD<< | |
1530 | "kink.="<<kink<< | |
1531 | "trackM.="<<trackM<< | |
1532 | "trackD.="<<trackD<< | |
1533 | "pm.="<<&pm<< //updated parameters | |
1534 | "pd.="<<&pd<< // updated parameters | |
1535 | "v0KF.="<<v0KF<< | |
1536 | "chi2="<<chi2<< | |
1537 | "\n"; | |
1538 | } | |
1539 | /* | |
1540 | TCut cutQ="chi2<10&&kink.fRr>90&&kink.fRr<220"; | |
1541 | TCut cutRD="20*sqrt(pd.fC[14])<abs(pd.fP[4])&&trackD.fTPCsignal>10&&trackD.fTPCsignalN>50"; | |
1542 | ||
1543 | */ | |
1544 | // | |
1545 | if (chi2>kChi2Cut) continue; | |
1546 | if (kink->GetR()<kMinR) continue; | |
1547 | if (kink->GetR()>kMaxR) continue; | |
1548 | if ((nclM+nclD)<kMinNcl) continue; | |
1549 | if (TMath::Abs(eta)>1) continue; | |
1550 | // | |
1551 | // | |
1552 | AliESDfriendTrack *friendTrackM = esdFriend->GetTrack(kink->GetIndex(0)); | |
1553 | AliESDfriendTrack *friendTrackD = esdFriend->GetTrack(kink->GetIndex(1)); | |
1554 | if (!friendTrackM) continue; | |
1555 | if (!friendTrackD) continue; | |
1556 | TObject *calibObject; | |
1557 | AliTPCseed *seedM = 0; | |
1558 | AliTPCseed *seedD = 0; | |
1559 | for (Int_t l=0;(calibObject=friendTrackM->GetCalibObject(l));++l) { | |
1560 | if ((seedM=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
1561 | } | |
1562 | for (Int_t l=0;(calibObject=friendTrackD->GetCalibObject(l));++l) { | |
1563 | if ((seedD=dynamic_cast<AliTPCseed*>(calibObject))) break; | |
1564 | } | |
1565 | } | |
1566 | } | |
1567 | ||
1568 | void AliTPCcalibGainMult::DumpHPT(const AliESDEvent * event){ | |
1569 | // | |
1570 | // Function to select the HPT tracks and events | |
1571 | // It is used to select event with HPT - list used later for the raw data downloading | |
1572 | // - and reconstruction | |
1573 | // Not actualy used for the calibration of the data | |
1574 | ||
1575 | TTreeSRedirector * pcstream = GetDebugStreamer(); | |
1576 | AliKFParticle::SetField(event->GetMagneticField()); | |
1577 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend")); | |
1578 | if (!esdFriend) { | |
1579 | //Printf("ERROR: esdFriend not available"); | |
1580 | return; | |
1581 | } | |
1582 | if (esdFriend->TestSkipBit()) return; | |
1583 | ||
1584 | Int_t ntracks=event->GetNumberOfTracks(); | |
1585 | // | |
1586 | for (Int_t i=0;i<ntracks;++i) { | |
1587 | // | |
1588 | AliESDtrack *track = event->GetTrack(i); | |
1589 | if (!track) continue; | |
1590 | if (track->Pt()<4) continue; | |
1591 | UInt_t status = track->GetStatus(); | |
1592 | // | |
1593 | AliExternalTrackParam * trackIn = (AliExternalTrackParam *)track->GetInnerParam(); | |
1594 | if (!trackIn) continue; | |
1595 | if ((status&AliESDtrack::kTPCrefit)==0) continue; | |
1596 | if ((status&AliESDtrack::kITSrefit)==0) continue; | |
1597 | AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i); | |
1598 | if (!friendTrack) continue; | |
1599 | AliExternalTrackParam * itsOut = (AliExternalTrackParam *)(friendTrack->GetITSOut()); | |
1600 | if (!itsOut) continue; | |
1601 | AliExternalTrackParam * itsOut2 = (AliExternalTrackParam *)(friendTrack->GetITSOut()->Clone()); | |
1602 | AliExternalTrackParam * tpcIn2 = (AliExternalTrackParam *)(trackIn->Clone()); | |
1603 | if (!itsOut2->Rotate(trackIn->GetAlpha())) continue; | |
1604 | //Double_t xmiddle=0.5*(itsOut2->GetX()+tpcIn2->GetX()); | |
1605 | Double_t xmiddle=(itsOut2->GetX()); | |
1606 | if (!itsOut2->PropagateTo(xmiddle,event->GetMagneticField())) continue; | |
1607 | if (!tpcIn2->PropagateTo(xmiddle,event->GetMagneticField())) continue; | |
1608 | // | |
1609 | AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam()); | |
1610 | if (!tpcInner) continue; | |
1611 | tpcInner->Rotate(track->GetAlpha()); | |
1612 | tpcInner->PropagateTo(track->GetX(),event->GetMagneticField()); | |
1613 | // | |
1614 | // tpc constrained | |
1615 | // | |
1616 | AliExternalTrackParam * tpcInnerC = (AliExternalTrackParam *)(track->GetTPCInnerParam()->Clone()); | |
1617 | if (!tpcInnerC) continue; | |
1618 | tpcInnerC->Rotate(track->GetAlpha()); | |
1619 | tpcInnerC->PropagateTo(track->GetX(),event->GetMagneticField()); | |
1620 | Double_t dz[2],cov[3]; | |
1621 | AliESDVertex *vtx= (AliESDVertex *)event->GetPrimaryVertex(); | |
1622 | ||
1623 | if (!tpcInnerC->PropagateToDCA(vtx, event->GetMagneticField(), 3, dz, cov)) continue; | |
1624 | Double_t covar[6]; vtx->GetCovMatrix(covar); | |
1625 | Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]}; | |
1626 | Double_t c[3]={covar[2],0.,covar[5]}; | |
1627 | // | |
1628 | Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c); | |
1629 | tpcInnerC->Update(p,c); | |
1630 | ||
1631 | if (pcstream){ | |
1632 | (*pcstream)<<"hpt"<< | |
1633 | "run="<<fRun<< | |
1634 | "time="<<fTime<< | |
1635 | "vertex="<<vtx<< | |
1636 | "bz="<<fMagF<< | |
1637 | "track.="<<track<< | |
1638 | "tpcInner.="<<tpcInner<< | |
1639 | "tpcInnerC.="<<tpcInnerC<< | |
1640 | "chi2C="<<chi2C<< | |
1641 | // | |
1642 | "its.="<<itsOut<< | |
1643 | "its2.="<<itsOut2<< | |
1644 | "tpc.="<<trackIn<< | |
1645 | "tpc2.="<<tpcIn2<< | |
1646 | "\n"; | |
1647 | } | |
1648 | } | |
1649 | } | |
1650 | ||
1651 | ||
1652 | ||
1653 | void AliTPCcalibGainMult::ProcessTOF(const AliESDEvent * event){ | |
1654 | // | |
1655 | // 1. Loop over tracks | |
1656 | // 2. Fit T0 | |
1657 | // 3. Sign positivelly identified tracks | |
1658 | // | |
1659 | const Double_t kMaxDelta=1000; | |
1660 | const Double_t kOrbit=50000; // distance in the time beween two orbits in the TOF time units - 50000=50 ns | |
1661 | const Double_t kMaxD=20; | |
1662 | const Double_t kRMS0=200; | |
1663 | const Double_t kMaxDCAZ=10; | |
1664 | AliESDVertex *vtx= (AliESDVertex *)event->GetPrimaryVertex(); | |
1665 | // | |
1666 | TTreeSRedirector * pcstream = GetDebugStreamer(); | |
1667 | AliKFParticle::SetField(event->GetMagneticField()); | |
1668 | AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend")); | |
1669 | if (!esdFriend) { | |
1670 | //Printf("ERROR: esdFriend not available"); | |
1671 | return; | |
1672 | } | |
1673 | //if (esdFriend->TestSkipBit()) return; | |
1674 | ||
1675 | Int_t ntracks=event->GetNumberOfTracks(); | |
1676 | // | |
1677 | Double_t deltaTPion[10000]; | |
1678 | Double_t medianT0=0; | |
1679 | Double_t meanT0=0; | |
1680 | Double_t rms=10000; | |
1681 | Int_t counter=0; | |
1682 | // | |
1683 | // Get Median time for pion hypothesy | |
1684 | // | |
1685 | for (Int_t iter=0; iter<3; iter++){ | |
1686 | counter=0; | |
1687 | for (Int_t i=0;i<ntracks;++i) { | |
1688 | // | |
1689 | AliESDtrack *track = event->GetTrack(i); | |
1690 | if (!track) continue; | |
1691 | if (!track->IsOn(AliESDtrack::kTIME)) continue; | |
1692 | if (TMath::Abs(track->GetZ())>kMaxDCAZ) continue; // remove overlaped events | |
1693 | if (TMath::Abs(track->GetTOFsignalDz())>kMaxD) continue; | |
1694 | Double_t times[1000]; | |
1695 | track->GetIntegratedTimes(times); | |
1696 | Int_t norbit=TMath::Nint((track->GetTOFsignal()-times[2])/kOrbit); | |
1697 | Double_t torbit=norbit*kOrbit; | |
1698 | if (iter==1 &&TMath::Abs(times[2]-times[3])<3*rms) continue; // skip umbigous points - kaon pion | |
1699 | // | |
1700 | Int_t indexBest=2; | |
1701 | if (iter>1){ | |
1702 | for (Int_t j=3; j<5; j++) | |
1703 | if (TMath::Abs(track->GetTOFsignal()-times[j]-torbit-medianT0)<TMath::Abs(track->GetTOFsignal()-times[j]-torbit-medianT0)) indexBest=j; | |
1704 | } | |
1705 | // | |
1706 | if (iter>0) if (TMath::Abs(track->GetTOFsignal()-times[indexBest]-torbit-medianT0)>3*(kRMS0+rms)) continue; | |
1707 | if (iter>0) if (TMath::Abs(track->GetTOFsignal()-times[indexBest]-torbit-medianT0)>kMaxDelta) continue; | |
1708 | deltaTPion[counter]=track->GetTOFsignal()-times[indexBest]-torbit; | |
1709 | counter++; | |
1710 | } | |
1711 | if (counter<2) return; | |
1712 | medianT0=TMath::Median(counter,deltaTPion); | |
1713 | meanT0=TMath::Median(counter,deltaTPion); | |
1714 | rms=TMath::RMS(counter,deltaTPion); | |
1715 | } | |
1716 | if (counter<3) return; | |
1717 | // | |
1718 | // Dump | |
1719 | // | |
1720 | for (Int_t i=0;i<ntracks;++i) { | |
1721 | // | |
1722 | AliESDtrack *track = event->GetTrack(i); | |
1723 | if (!track) continue; | |
1724 | if (!track->IsOn(AliESDtrack::kTIME)) continue; | |
1725 | if (TMath::Abs(track->GetZ())>kMaxDCAZ) continue; //remove overlapped events | |
1726 | if (TMath::Abs(track->GetTOFsignalDz())>kMaxD) continue; | |
1727 | Double_t times[1000]; | |
1728 | track->GetIntegratedTimes(times); | |
1729 | Int_t norbit=TMath::Nint((track->GetTOFsignal()-times[2])/kOrbit); | |
1730 | Double_t torbit=norbit*kOrbit; | |
1731 | if (rms<=0) continue; | |
1732 | // | |
1733 | Double_t tPion = (track->GetTOFsignal()-times[2]-medianT0-torbit); | |
1734 | Double_t tKaon = (track->GetTOFsignal()-times[3]-medianT0-torbit); | |
1735 | Double_t tProton= (track->GetTOFsignal()-times[4]-medianT0-torbit); | |
1736 | Double_t tElectron= (track->GetTOFsignal()-times[0]-medianT0-torbit); | |
1737 | // | |
1738 | Bool_t isPion = (TMath::Abs(tPion/rms)<6) && TMath::Abs(tPion)<(TMath::Min(TMath::Abs(tKaon), TMath::Abs(tProton))-rms); | |
1739 | Bool_t isKaon = (TMath::Abs(tKaon/rms)<3) && TMath::Abs(tKaon)<0.2*(TMath::Min(TMath::Abs(tPion), TMath::Abs(tProton))-3*rms); | |
1740 | Bool_t isProton = (TMath::Abs(tProton/rms)<6) && TMath::Abs(tProton)<0.5*(TMath::Min(TMath::Abs(tKaon), TMath::Abs(tPion))-rms); | |
1741 | Bool_t isElectron = (TMath::Abs(tElectron/rms)<6) && TMath::Abs(tElectron)<0.2*(TMath::Min(TMath::Abs(tKaon), TMath::Abs(tPion))-rms) &&TMath::Abs(tElectron)<0.5*(TMath::Abs(tPion)-rms); | |
1742 | ||
1743 | if (isPion) (*fPIDMatrix)(i,2)+=1; | |
1744 | if (isKaon) (*fPIDMatrix)(i,3)+=1; | |
1745 | if (isProton) (*fPIDMatrix)(i,4)+=1; | |
1746 | // if (isElectron) (*fPIDMatrix)(i,0)+=1; | |
1747 | // | |
1748 | if (pcstream){ | |
1749 | // debug streamer to dump the information | |
1750 | (*pcstream)<<"tof"<< | |
1751 | "isPion="<<isPion<< | |
1752 | "isKaon="<<isKaon<< | |
1753 | "isProton="<<isProton<< | |
1754 | "isElectron="<<isElectron<< | |
1755 | // | |
1756 | "counter="<<counter<< | |
1757 | "torbit="<<torbit<< | |
1758 | "norbit="<<norbit<< | |
1759 | "medianT0="<<medianT0<< | |
1760 | "meanT0="<<meanT0<< | |
1761 | "rmsT0="<<rms<< | |
1762 | "track.="<<track<< | |
1763 | "vtx.="<<vtx<< | |
1764 | "\n"; | |
1765 | } | |
1766 | ||
1767 | } | |
1768 | /* | |
1769 | tof->SetAlias("isProton","(abs(track.fTOFsignal-track.fTrackTime[4]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)-rmsT0))"); | |
1770 | tof->SetAlias("isPion","(abs(track.fTOFsignal-track.fTrackTime[2]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)-rmsT0))"); | |
1771 | tof->SetAlias("isKaon","(abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[2]-medianT0-torbit)-rmsT0))&&(abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[4]-medianT0-torbit)-rmsT0))"); | |
1772 | ||
1773 | */ | |
1774 | ||
1775 | } | |
1776 | ||
1777 | ||
1778 | // void AliTPCcalibGainMult::Terminate(){ | |
1779 | // // | |
1780 | // // Terminate function | |
1781 | // // call base terminate + Eval of fitters | |
1782 | // // | |
1783 | // Info("AliTPCcalibGainMult","Terminate"); | |
1784 | // TTreeSRedirector *pcstream = GetDebugStreamer(); | |
1785 | // if (pcstream){ | |
1786 | // TTreeStream &stream = (*pcstream)<<"dump"; | |
1787 | // TTree* tree = stream.GetTree(); | |
1788 | // if (tree) if ( tree->GetEntries()>0){ | |
1789 | // TObjArray *array = tree->GetListOfBranches(); | |
1790 | // for (Int_t i=0; i<array->GetEntries(); i++) {TBranch * br = (TBranch *)array->At(i); br->SetAddress(0);} | |
1791 | // gDirectory=gROOT; | |
1792 | // fdEdxTree=tree->CloneTree(10000); | |
1793 | // } | |
1794 | // } | |
1795 | // AliTPCcalibBase::Terminate(); | |
1796 | // } | |
1797 |