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