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