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