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