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