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