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