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