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