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