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