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[6]    = { 400, 400,    4,   20,   50, 100};
 151   Double_t xminPadEqual[6] = { 0.0, 0.0, -0.5,    0, -250,   0};
 152   Double_t xmaxPadEqual[6] = { 2.0, 2.0,  3.5, 13000,  250,   3};
 153   TString axisNamePadEqual[6]   = {"dEdxRatioMax","dEdxRatioTot","padType","mult","driftlength", "1_pt"};
 154   TString axisTitlePadEqual[6]  = {"dEdx_padRegion/mean_dEdx Qmax", "dEdx_padRegion/mean_dEdx Qtot","padType","mult","driftlength", "1/pt"};
 155   //
 156   fHistPadEqual = new THnSparseF("fHistPadEqual","0:dEdx_pad/dEdx_mean, 1:pad, 2:mult, 3:drift, 4:1/pt", 6, binsPadEqual, xminPadEqual, xmaxPadEqual);
 157   for (Int_t iaxis=0; iaxis<6;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[6] = {signalArrayMax[ipad]/meanMax, signalArrayTot[ipad]/meanTot, ipad, nContributors, meanDrift, track->OneOverPt()};
 377         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   TIterator* iter = li->MakeIterator();
 466   AliTPCcalibGainMult* cal = 0;
 467
 468   while ((cal = (AliTPCcalibGainMult*)iter->Next())) {
 469     if (!cal->InheritsFrom(AliTPCcalibGainMult::Class())) {
 470       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 471       return -1;
 472     }
 473
 474     if (cal->GetHistNTracks()) fHistNTracks->Add(cal->GetHistNTracks());
 475     if (cal->GetHistClusterShape()) fHistClusterShape->Add(cal->GetHistClusterShape());
 476     if (cal->GetHistQA()) fHistQA->Add(cal->GetHistQA());
 477     if (cal->GetHistGainSector()) fHistGainSector->Add(cal->GetHistGainSector());
 478     if (cal->GetHistPadEqual()) fHistPadEqual->Add(cal->GetHistPadEqual());
 479     if (cal->GetHistGainMult()) fHistGainMult->Add(cal->GetHistGainMult());
 480     if (cal->fHistdEdxMap){
 481       if (fHistdEdxMap) fHistdEdxMap->Add(cal->fHistdEdxMap);
 482     }
 483     if (cal->fHistdEdxMax){
 484       if (fHistdEdxMax) fHistdEdxMax->Add(cal->fHistdEdxMax);
 485     }
 486     if (cal->fHistdEdxTot){
 487       if (fHistdEdxTot) fHistdEdxTot->Add(cal->fHistdEdxTot);
 488     }
 489     //
 490     // Originally we tireied to write the tree to the same file as other calibration components
 491     // We failed in merging => therefore this optio  was disabled
 492     //
 493     //    if (cal->fdEdxTree && cal->fdEdxTree->GetEntries()>0) {
 494     //       if (fdEdxTree) {
 495     //  const Int_t kMax=100000;
 496     //  Int_t entriesSum = (Int_t)fdEdxTree->GetEntries();
 497     //  Int_t entriesCurrent = (Int_t)cal->fdEdxTree->GetEntries();
 498     //  Int_t entriesCp=TMath::Min((Int_t) entriesCurrent*(kMax*entriesSum),entriesCurrent);
 499     // //       cal->fdEdxTree->SetBranchStatus("track*",0);
 500     // //       cal->fdEdxTree->SetBranchStatus("vertex*",0);
 501     // //       cal->fdEdxTree->SetBranchStatus("tpcOut*",0);
 502     // //       cal->fdEdxTree->SetBranchStatus("vec*",0);
 503     // //       fdEdxTree->SetBranchStatus("track*",0);
 504     // //       fdEdxTree->SetBranchStatus("vertex*",0);
 505     // //       fdEdxTree->SetBranchStatus("tpcOut*",0);
 506     // //       fdEdxTree->SetBranchStatus("vec*",0);
 507     //  fdEdxTree->Print();
 508     //  fdEdxTree->Dump();
 509     //  fdEdxTree->GetEntry(0);
 510     //  for (Int_t i=0; i<entriesCurrent; i++){
 511     //    cal->fdEdxTree->CopyAddresses(fdEdxTree);
 512     //    cal->fdEdxTree->GetEntry(i);
 513     //    fdEdxTree->Fill();
 514     //  }
 515     //  TObjArray *brarray =  cal->fdEdxTree->GetListOfBranches();
 516     //  for (Int_t i=0; i<brarray->GetEntries(); i++) {TBranch * br = (TBranch *)brarray->At(i); br->SetAddress(0); }
 517     //       }
 518     //       else{
 519     //  fdEdxTree = (TTree*)(cal->fdEdxTree->Clone());
 520     //  TObjArray *brarray =  fdEdxTree->GetListOfBranches();
 521     //  for (Int_t i=0; i<brarray->GetEntries(); i++) {TBranch * br = (TBranch *)brarray->At(i); br->SetAddress(0);}
 522     //       }
 523     //}
 524
 525   }
 526
 527   return 0;
 528
 529 }
 530
 531
 532
 533
 534
 535 void AliTPCcalibGainMult::UpdateGainMap() {
 536   //
 537   // read in the old gain map and scale it appropriately...
 538   //
 539   /*
 540   gSystem->Load("libANALYSIS");
 541   gSystem->Load("libTPCcalib");
 542   //
 543   TFile jj("Run0_999999999_v1_s0.root");
 544   AliTPCCalPad * pad = AliCDBEntry->GetObject()->Clone();
 545   TFile hh("output.root");
 546   AliTPCcalibGainMult * gain = calibTracksGain;
 547   TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,1);
 548   //
 549   TObjArray arr;
 550   histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
 551   TH1D * meanGainSec = arr->At(1);
 552   Double_t gainsIROC[36];
 553   Double_t gainsOROC[36];
 554   Double_t gains[72];
 555   //
 556   for(Int_t isec = 1; isec < meanGainSec->GetNbinsX() + 1; isec++) {
 557     cout << isec << " " << meanGainSec->GetXaxis()->GetBinCenter(isec) << " " <<meanGainSec->GetBinContent(isec) << endl;
 558     gains[isec-1] = meanGainSec->GetBinContent(isec);
 559     if (isec < 37) {
 560       gainsIROC[isec-1] = meanGainSec->GetBinContent(isec);
 561     } else {
 562       gainsOROC[isec - 36 -1] = meanGainSec->GetBinContent(isec);
 563     }
 564   }
 565   Double_t meanIroc = TMath::Mean(36, gainsIROC);
 566   Double_t meanOroc = TMath::Mean(36, gainsIROC);
 567   for(Int_t i = 0; i < 36; i++) gains[i] /= meanIroc;
 568   for(Int_t i = 36; i < 72; i++) gains[i] /= meanOroc;
 569   //
 570   for(Int_t i = 0; i < 72; i++) {
 571     AliTPCCalROC * chamber = pad->GetCalROC(i);
 572     chamber->Multiply(gains[i]);
 573     cout << i << " "<< chamber->GetMean() << endl;
 574   }
 575   //
 576   // update the OCDB
 577   //
 578   AliCDBMetaData *metaData= new AliCDBMetaData();
 579   metaData->SetObjectClassName("AliTPCCalPad");
 580   metaData->SetResponsible("Alexander Kalweit");
 581   metaData->SetBeamPeriod(1);
 582   metaData->SetAliRootVersion("04-19-05"); //root version
 583   metaData->SetComment("New gain map for 1600V OROC gain increase and equalization. Valid for runs starting after technical stop beginning of September.");
 584   AliCDBId id1("TPC/Calib/GainFactorDedx", 131541, AliCDBRunRange::Infinity()); // important: new gain runs here..
 585   AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///d/alice05/akalweit/projects/OCDBupdate/HighGain_2010-09-03/OCDB/");
 586   gStorage->Put(pad, id1, metaData);
 587   */
 588
 589 }
 590
 591 void AliTPCcalibGainMult::UpdateClusterParam() {
 592   //
 593   //
 594   //
 595   /*
 596   gSystem->Load("libANALYSIS");
 597   gSystem->Load("libTPCcalib");
 598   //
 599   TFile ff("OldClsParam.root");
 600   AliTPCClusterParam * param = AliCDBEntry->GetObject()->Clone();
 601
 602   TFile hh("output.root");
 603   AliTPCcalibGainMult * gain = calibGainMult;
 604   TH2D * histGainSec = gain->GetHistGainSector()->Projection(0,2);
 605   TObjArray arr;
 606   histGainSec->FitSlicesY(0, 0, -1, 0, "QNR", &arr);
 607   histGainSec->Draw("colz");
 608   TH1D * fitVal = arr.At(1);
 609   fitVal->Draw("same");
 610   param->GetQnormCorrMatrix()->Print();
 611   param->GetQnormCorrMatrix()(0,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qtot
 612   param->GetQnormCorrMatrix()(1,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qtot
 613   param->GetQnormCorrMatrix()(2,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qtot
 614   //
 615   param->GetQnormCorrMatrix()(3,5) *= fitVal->GetBinContent(1)/fitVal->GetBinContent(1); // short pads Qmax -> scaling assumed
 616   param->GetQnormCorrMatrix()(4,5) *= fitVal->GetBinContent(2)/fitVal->GetBinContent(1); // med pads Qmax -> scaling assumed
 617   param->GetQnormCorrMatrix()(5,5) *= fitVal->GetBinContent(3)/fitVal->GetBinContent(1); // long pads Qmax -> scaling assumed
 618   //
 619   TFile jj("newClusterParam.root","RECREATE");
 620   param->Write();
 621   param->GetQnormCorrMatrix()->Print();
 622   //
 623   // update the OCDB
 624   //
 625   AliCDBMetaData *metaData= new AliCDBMetaData();
 626   metaData->SetObjectClassName("AliTPCClusterParam");
 627   metaData->SetResponsible("Alexander Kalweit");
 628   metaData->SetBeamPeriod(1);
 629   metaData->SetAliRootVersion("04-19-04"); //root version
 630   metaData->SetComment("1600V OROC / hard thres. / new algorithm");
 631   AliCDBId id1("TPC/Calib/ClusterParam", 0, AliCDBRunRange::Infinity()); // important: new gain runs here..
 632   AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage("local:///lustre/alice/akalweit/baseline/CalibrationEntries/OldThres_NewAlgo_PP");
 633   gStorage->Put(param, id1, metaData);
 634   */
 635
 636
 637 }
 638
 639
 640 void AliTPCcalibGainMult::DumpTrack(AliESDtrack * track, AliESDfriendTrack *ftrack, AliTPCseed * seed, Int_t index){
 641   //
 642   // dump interesting tracks
 643   // 1. track at MIP region
 644   // 2. highly ionizing protons
 645   // pidType: 0 - unselected
 646   //          1 - TOF
 647   //          2 - V0
 648   //          4 - Cosmic
 649   //          or of value
 650   //
 651   const Int_t    kMax=10000;
 652   const Int_t    kMinRows=80;
 653   const Double_t kDCAcut=30;
 654   //
 655   // Bethe Bloch paramerization
 656   //
 657   Double_t kp1= 0.0851148;
 658   Double_t kp2= 9.25771;
 659   Double_t kp3= 2.6558e-05;
 660   Double_t kp4= 2.32742;
 661   Double_t kp5= 1.83039;
 662   if (fBBParam){
 663     kp1=(*fBBParam)[0];
 664     kp2=(*fBBParam)[1];
 665     kp3=(*fBBParam)[2];
 666     kp4=(*fBBParam)[3];
 667     kp5=(*fBBParam)[4];
 668   }
 669   //
 670   static const AliTPCROC *roc = AliTPCROC::Instance();
 671   static const TDatabasePDG *pdg = TDatabasePDG::Instance();
 672
 673   Int_t nclITS   = track->GetNcls(0);
 674   Int_t ncl   = track->GetTPCncls();
 675   Double_t ncl21 = track->GetTPCClusterInfo(3,1);
 676   Double_t ncl20 = track->GetTPCClusterInfo(3,0);
 677   //
 678   if (!seed) return;
 679   if (ncl21<kMinRows) return;
 680   static Int_t counter=0;
 681   static Int_t counterHPT=0;
 682   //
 683   static TH1F     *hisBB=new TH1F("hisBB","hisBB",20,0.1,1.00);  // bethe bloch histogram  =
 684   //                                                                 used to cover more homegenously differnt dEdx regions
 685   static Double_t massPi = pdg->GetParticle("pi-")->Mass();      //
 686   static Double_t massK  = pdg->GetParticle("K-")->Mass();
 687   static Double_t massP  = pdg->GetParticle("proton")->Mass();
 688   static Double_t massE  = pdg->GetParticle("e-")->Mass();
 689   static Double_t massMuon  = pdg->GetParticle("mu-")->Mass();
 690   static Double_t radius0= roc->GetPadRowRadiiLow(roc->GetNRows(0)/2);
 691   static Double_t radius1= roc->GetPadRowRadiiUp(30);
 692   static Double_t radius2= roc->GetPadRowRadiiUp(roc->GetNRows(36)-15);
 693
 694   AliESDVertex *vertex= (AliESDVertex *)fCurrentEvent->GetPrimaryVertex();
 695   //
 696   // Estimate current MIP position -
 697   //
 698   static Double_t mipArray[kMax];               // mip array
 699   static Int_t    counterMIP0=0;
 700   static Double_t    medianMIP0=100000;         // current MIP median position - estimated after some mimnimum number of MIP tracks
 701
 702   if (TMath::Abs(track->GetP()-0.5)<0.1&&track->GetTPCsignal()/medianMIP0<1.5){
 703     mipArray[counterMIP0%kMax]= track->GetTPCsignal();
 704     counterMIP0++;
 705     if (counterMIP0>10) medianMIP0=TMath::Median(counterMIP0%kMax, mipArray);
 706   }
 707   // the PID as defiend from the external sources
 708   //
 709   Int_t isElectron   =  TMath::Nint((*fPIDMatrix)(index,0));
 710   Int_t isMuon       =  TMath::Nint((*fPIDMatrix)(index,1));
 711   Int_t isPion       =  TMath::Nint((*fPIDMatrix)(index,2));
 712   Int_t isKaon       =  TMath::Nint((*fPIDMatrix)(index,3));
 713   Int_t isProton     =  TMath::Nint((*fPIDMatrix)(index,4));
 714   Float_t dca[2];
 715   track->GetImpactParameters(dca[0],dca[1]);
 716   //
 717   if ( (isMuon==0 && isElectron==0)  && (TMath::Sqrt(dca[0]*dca[0]+dca[1]*dca[1])>kDCAcut) ) return;
 718   Double_t normdEdx= track->GetTPCsignal()/(medianMIP0); // TPC signal normalized to the MIP
 719   //
 720   AliExternalTrackParam * trackIn  = (AliExternalTrackParam *)track->GetInnerParam();
 721   AliExternalTrackParam * trackOut = (AliExternalTrackParam *)track->GetOuterParam();
 722   AliExternalTrackParam * tpcOut   = (AliExternalTrackParam *)ftrack->GetTPCOut();
 723   if (!trackIn) return;
 724   if (!trackOut) return;
 725   if (!tpcOut) return;
 726   if (trackIn->GetZ()*trackOut->GetZ()<0) return;  // remove crossing tracks
 727   //
 728   // calculate local and global angle
 729   Int_t side = (trackIn->GetZ()>0)? 1:-1;
 730   Double_t tgl=trackIn->GetTgl();
 731   Double_t gangle[3]={0,0,0};
 732   Double_t langle[3]={0,0,0};
 733   Double_t length[3]={0,0,0};
 734   Double_t pxpypz[3]={0,0,0};
 735   Double_t bz=fMagF;
 736   trackIn->GetXYZAt(radius0,bz,pxpypz);            // get the global position  at the middle of the IROC
 737   gangle[0]=TMath::ATan2(pxpypz[1],pxpypz[0]);     // global angle IROC
 738   trackIn->GetXYZAt(radius1,bz,pxpypz);            // get the global position at the middle of the OROC - medium pads
 739   gangle[1]=TMath::ATan2(pxpypz[1],pxpypz[0]);     // global angle OROC
 740   trackOut->GetXYZAt(radius2,bz,pxpypz);           // get the global position at the middle of OROC - long pads
 741   gangle[2]=TMath::ATan2(pxpypz[1],pxpypz[0]);
 742   //
 743   trackIn->GetPxPyPzAt(radius0,bz,pxpypz);               //get momentum vector
 744   langle[0]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[0];  //local angle between padrow and track IROC
 745   trackIn->GetPxPyPzAt(radius1,bz,pxpypz);
 746   langle[1]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[1];
 747   trackOut->GetPxPyPzAt(radius2,bz,pxpypz);               //                                     OROC medium
 748   langle[2]=TMath::ATan2(pxpypz[1],pxpypz[0])-gangle[2];
 749   for (Int_t i=0; i<3; i++){
 750     if (langle[i]>TMath::Pi())  langle[i]-=TMath::TwoPi();
 751     if (langle[i]<-TMath::Pi()) langle[i]+=TMath::TwoPi();
 752     length[i]=TMath::Sqrt(1+langle[i]*langle[i]+tgl*tgl);  // the tracklet length
 753   }
 754   //
 755   // Select the kaons and Protons which are "isolated" in TPC dedx curve
 756   //
 757   //
 758   Double_t dedxP = AliExternalTrackParam::BetheBlochAleph(track->GetInnerParam()->GetP()/massP,kp1,kp2,kp3,kp4,kp5);
 759   Double_t dedxK = AliExternalTrackParam::BetheBlochAleph(track->GetInnerParam()->GetP()/massK,kp1,kp2,kp3,kp4,kp5);
 760   if (dedxP>2 || dedxK>2){
 761     if (track->GetP()<1.2 && normdEdx>1.8&&counterMIP0>10){ // not enough from TOF and V0s triggered by high dedx
 762       // signing the Proton  and kaon - using the "bitmask" bit 1 and 2 is dedicated for V0s and TOF selected
 763       if ( TMath::Abs(normdEdx/dedxP-1)<0.3)  isProton+=4;
 764       if ( TMath::Abs(normdEdx/dedxK-1)<0.3)  isKaon+=4;
 765       if (normdEdx/dedxK>1.3) isProton+=8;
 766       if (normdEdx/dedxP<0.7) isKaon+=8;
 767     }
 768   }
 769   //
 770   //
 771   //
 772   Double_t mass = 0;
 773   Bool_t isHighPt = ((TMath::Power(1/track->Pt(),4)*gRandom->Rndm())<0.005);  // rnadomly selected HPT tracks
 774   // there are selected for the QA of the obtained calibration
 775   Bool_t isMIP    =  TMath::Abs(track->GetInnerParam()->P()-0.4)<0.005&&(counter<kMax); //
 776   // REMINDER - it is not exactly MIP - we select the regtion where the Kaon and Electrons are well separated
 777
 778   if (isElectron>0) mass = massE;
 779   if (isProton>0)   mass = massP;
 780   if (isKaon>0)     mass = massK;
 781   if (isMuon>0)     mass = massMuon;
 782   if (isPion>0)     mass = massPi;
 783   if (isHighPt)     mass = massPi;  //assign mass of pions
 784   if (isMIP&&track->GetTPCsignal()/medianMIP0<1.5)   mass = massPi;  //assign mass of pions
 785   if (mass==0)      return;
 786   //
 787   // calculate expected dEdx
 788   Double_t dedxDef= 0;
 789   Double_t dedxDefPion= 0,dedxDefProton=0, dedxDefKaon=0;
 790   Double_t pin=trackIn->GetP();
 791   Double_t pout=trackOut->GetP();
 792   Double_t p=(pin+pout)*0.5;  // momenta as the mean between tpc momenta at inner and outer wall of the TPC
 793   if (mass>0) dedxDef = AliExternalTrackParam::BetheBlochAleph(p/mass,kp1,kp2,kp3,kp4,kp5);
 794   dedxDefPion = AliExternalTrackParam::BetheBlochAleph(p/massPi,kp1,kp2,kp3,kp4,kp5);
 795   dedxDefProton = AliExternalTrackParam::BetheBlochAleph(p/massP,kp1,kp2,kp3,kp4,kp5);
 796   dedxDefKaon = AliExternalTrackParam::BetheBlochAleph(p/massK,kp1,kp2,kp3,kp4,kp5);
 797   //
 798   // dEdx Truncated mean vectros with differnt tuncation
 799   // 11 truncations array -  0-10  - 0~50%  11=100%
 800   // 3 Regions            -  IROC,OROC0, OROC1
 801   // 2 Q                  -  total charge and max charge
 802   // Log                  -  Logarithmic mean used
 803   // Up/Dwon              -  Upper half or lower half of truncation used
 804   // RMS                  -  rms of the distribction (otherwise truncated mean)
 805   // M2 suffix            -  second moment ( truncated)
 806   TVectorF truncUp(11);
 807   TVectorF truncDown(11);
 808   TVectorF vecAllMax(11);
 809   TVectorF vecIROCMax(11);
 810   TVectorF vecOROCMax(11);
 811   TVectorF vecOROC0Max(11);
 812   TVectorF vecOROC1Max(11);
 813   //
 814   TVectorF vecAllTot(11);
 815   TVectorF vecIROCTot(11);
 816   TVectorF vecOROCTot(11);
 817   TVectorF vecOROC0Tot(11);
 818   TVectorF vecOROC1Tot(11);
 819   //
 820   TVectorF vecAllTotLog(11);
 821   TVectorF vecIROCTotLog(11);
 822   TVectorF vecOROCTotLog(11);
 823   TVectorF vecOROC0TotLog(11);
 824   TVectorF vecOROC1TotLog(11);
 825   //
 826   TVectorF vecAllTotUp(11);
 827   TVectorF vecIROCTotUp(11);
 828   TVectorF vecOROCTotUp(11);
 829   TVectorF vecOROC0TotUp(11);
 830   TVectorF vecOROC1TotUp(11);
 831   //
 832   TVectorF vecAllTotDown(11);
 833   TVectorF vecIROCTotDown(11);
 834   TVectorF vecOROCTotDown(11);
 835   TVectorF vecOROC0TotDown(11);
 836   TVectorF vecOROC1TotDown(11);
 837
 838   TVectorF vecAllTotRMS(11);
 839   TVectorF vecIROCTotRMS(11);
 840   TVectorF vecOROCTotRMS(11);
 841   TVectorF vecOROC0TotRMS(11);
 842   TVectorF vecOROC1TotRMS(11);
 843   //
 844   TVectorF vecAllTotM2(11);
 845   TVectorF vecIROCTotM2(11);
 846   TVectorF vecOROCTotM2(11);
 847   TVectorF vecOROC0TotM2(11);
 848   TVectorF vecOROC1TotM2(11);
 849   //
 850   TVectorF vecAllTotMS(11);
 851   TVectorF vecIROCTotMS(11);
 852   TVectorF vecOROCTotMS(11);
 853   TVectorF vecOROC0TotMS(11);
 854   TVectorF vecOROC1TotMS(11);
 855   //
 856   // Differnt number of clusters definitions - in separate regions of the TPC
 857   // 20  -    ratio - found/findabel
 858   // 21  -    number of clusters used for given dEdx calculation
 859   //
 860   // suffix - 3 or 4 -  number of padrows before and after given row to define findable row
 861   //
 862   Double_t ncl20All  = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3);
 863   Double_t ncl20IROC = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3);
 864   Double_t ncl20OROC = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3);
 865   Double_t ncl20OROC0= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3);
 866   Double_t ncl20OROC1= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3);
 867   //
 868   Double_t ncl20All4  = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3,4);
 869   Double_t ncl20IROC4 = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3,4);
 870   Double_t ncl20OROC4 = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3,4);
 871   Double_t ncl20OROC04= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3,4);
 872   Double_t ncl20OROC14= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3,4);
 873   //
 874   Double_t ncl20All3  = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,3,3);
 875   Double_t ncl20IROC3 = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,3,3);
 876   Double_t ncl20OROC3 = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,3,3);
 877   Double_t ncl20OROC03= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,3,3);
 878   Double_t ncl20OROC13= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,3,3);
 879   //
 880   Double_t ncl21All  = seed->CookdEdxAnalytical(0.0,1, 1 ,0,159,2);
 881   Double_t ncl21IROC = seed->CookdEdxAnalytical(0.,1, 1 ,0,63,2);
 882   Double_t ncl21OROC = seed->CookdEdxAnalytical(0.,1, 1 ,64,159,2);
 883   Double_t ncl21OROC0= seed->CookdEdxAnalytical(0.,1, 1 ,64,128,2);
 884   Double_t ncl21OROC1= seed->CookdEdxAnalytical(0.,1, 1 ,129,159,2);
 885   // calculate truncated dEdx - mean rms M2 ...
 886   Int_t ifrac=0;
 887   for (Int_t ifracDown=0; ifracDown<1; ifracDown++){
 888     for (Int_t ifracUp=0; ifracUp<11; ifracUp++){
 889       Double_t fracDown = 0.0+Double_t(ifracDown)*0.05;
 890       Double_t fracUp = 0.5+Double_t(ifracUp)*0.05;
 891       vecAllMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,0,159,0);
 892       vecIROCMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,0,63,0);
 893       vecOROCMax[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,64,159,0);
 894       vecOROC0Max[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,64,128,0);
 895       vecOROC1Max[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 1 ,129,159,0);
 896       //
 897       vecAllTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,0);
 898       vecIROCTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,0);
 899       vecOROCTot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,0);
 900       vecOROC0Tot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,0);
 901       vecOROC1Tot[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,0);
 902       //
 903       vecAllTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,0,2,1);
 904       vecIROCTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,0,2,1);
 905       vecOROCTotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,0,2,1);
 906       vecOROC0TotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,0,2,1);
 907       vecOROC1TotLog[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,0,2,1);
 908       //
 909       vecAllTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,4,2,1);
 910       vecIROCTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,4,2,1);
 911       vecOROCTotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,4,2,1);
 912       vecOROC0TotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,4,2,1);
 913       vecOROC1TotUp[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,4,2,1);
 914       //
 915       vecAllTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,5,2,1);
 916       vecIROCTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,5,2,1);
 917       vecOROCTotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,5,2,1);
 918       vecOROC0TotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,5,2,1);
 919       vecOROC1TotDown[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,5,2,1);
 920       //
 921       vecAllTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,1,2,0);
 922       vecIROCTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,1,2,0);
 923       vecOROCTotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,1,2,0);
 924       vecOROC0TotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,1,2,0);
 925       vecOROC1TotRMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,1,2,0);
 926       //
 927       vecAllTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,6,2,1);
 928       vecIROCTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,6,2,1);
 929       vecOROCTotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,6,2,1);
 930       vecOROC0TotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,6,2,1);
 931       vecOROC1TotM2[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,6,2,1);
 932       //
 933       vecAllTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,159,8,2,1);
 934       vecIROCTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,0,63,8,2,1);
 935       vecOROCTotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,159,8,2,1);
 936       vecOROC0TotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,64,128,8,2,1);
 937       vecOROC1TotMS[ifrac]= seed->CookdEdxAnalytical(fracDown,fracUp, 0 ,129,159,8,2,1);
 938       truncUp[ifrac]=fracUp;
 939       truncDown[ifrac]=fracDown;
 940       ifrac++;
 941     }
 942   }
 943   //
 944   // Fill histograms
 945   //
 946   if ((isKaon||isProton||isPion||isElectron||isMIP||isMuon&&(!isHighPt)) && dedxDef>0) {
 947     //
 948     Int_t ncont = vertex->GetNContributors();
 949     for (Int_t ipad=0; ipad<3; ipad++){
 950       // histogram with enahanced phi granularity - to make gain phi maps
 951       Double_t xxx[4]={0,gangle[ipad],1./dedxDef,ipad*2+((side>0)?0:1)};
 952       Double_t nclR=0;
 953       if (ipad==0)  {xxx[0]=vecIROCTot[4]/medianMIP0; nclR=ncl21IROC/63.;}
 954       if (ipad==1)  {xxx[0]=vecOROC0Tot[4]/medianMIP0;nclR=ncl21OROC0/63.;}
 955       if (ipad==2)  {xxx[0]=vecOROC1Tot[4]/medianMIP0;nclR=ncl21OROC1/32.;}
 956       xxx[0]/=dedxDef;
 957       if (xxx[0]>0) xxx[0]=1./xxx[0];
 958       if (TMath::Abs(langle[ipad])<0.25&&nclR>0.4)  fHistdEdxMap->Fill(xxx);
 959     }
 960     for (Int_t ipad=0; ipad<3; ipad++){
 961       //
 962       // this are histogram to define  overall main gain correction
 963       // Maybe dead end - we can not put all info which can be used into the THnSparse
 964       // It is keeped there for educational point of view
 965       //
 966       Double_t xxx[6]={0,1./dedxDef, TMath::Abs(langle[ipad]), TMath::Abs(tgl), ncont, ipad };
 967       if (ipad==0)  {xxx[0]=vecIROCTot[4]/medianMIP0;}
 968       if (ipad==1)  {xxx[0]=vecOROC0Tot[4]/medianMIP0;}
 969       if (ipad==2)  {xxx[0]=vecOROC1Tot[4]/medianMIP0;}
 970       xxx[0]/=dedxDef;
 971       if (xxx[0]>0) xxx[0]=1./xxx[0];
 972       if (xxx[0]>0) fHistdEdxTot->Fill(xxx);
 973       if (ipad==0)  {xxx[0]=vecIROCMax[4]/medianMIP0;}
 974       if (ipad==1)  {xxx[0]=vecOROC0Max[4]/medianMIP0;}
 975       if (ipad==2)  {xxx[0]=vecOROC1Max[4]/medianMIP0;}
 976       xxx[0]=dedxDef;
 977       if (xxx[0]>0) xxx[0]=1./xxx[0];
 978       fHistdEdxMax->Fill(xxx);
 979     }
 980   }
 981   //
 982   // Downscale  selected tracks before filling the tree
 983   //
 984   Bool_t isSelected = kFALSE;
 985   //
 986   if (isKaon||isProton||isPion||isElectron||isMIP||isMuon) isSelected=kTRUE;
 987   isHighPt = kFALSE;
 988   if (!isSelected) isHighPt = ((TMath::Power(1/track->Pt(),4)*gRandom->Rndm())<0.005);
 989   if (counter>kMax && ((1/track->Pt()*gRandom->Rndm())>kMax/counter)) return;
 990   isSelected|=isHighPt;
 991   //
 992   //
 993   //
 994   // Equalize statistic in BB bins - special care about pions
 995   Int_t entriesBB = (Int_t)hisBB->GetEntries();
 996   if ((isElectron==0 &&isMuon==0 && p<2.) && entriesBB>20 &&dedxDef>0.01){
 997     Int_t bin = hisBB->GetXaxis()->FindBin(1./dedxDef);
 998     Double_t cont = hisBB->GetBinContent(bin);
 999     Double_t mean =(entriesBB)/20.;
1000     if ((isPion>0)  && gRandom->Rndm()*cont > 0.1*mean) return;
1001     if ((isPion==0) && gRandom->Rndm()*cont > 0.25*mean) return;
1002   }
1003   if (!isSelected) return;
1004   if (dedxDef>0.01) hisBB->Fill(1./dedxDef);
1005   //
1006   if (isHighPt) counterHPT++;
1007   counter++;
1008   //
1009   TTreeSRedirector * pcstream =  GetDebugStreamer();
1010   Double_t ptrel0 = AliTPCcalibDB::GetPTRelative(fTime,fRun,0);
1011   Double_t ptrel1 = AliTPCcalibDB::GetPTRelative(fTime,fRun,1);
1012   Int_t sectorIn   = Int_t(18+9*(trackIn->GetAlpha()/TMath::Pi()))%18;
1013   Int_t sectorOut  = Int_t(18+9*(trackOut->GetAlpha()/TMath::Pi()))%18;
1014   //
1015   if (pcstream){
1016     (*pcstream)<<"dump"<<
1017       "vertex.="<<vertex<<
1018       "bz="<<fMagF<<
1019       "ptrel0="<<ptrel0<<
1020       "ptrel1="<<ptrel1<<
1021       "sectorIn="<<sectorIn<<
1022       "sectorOut="<<sectorOut<<
1023       "side="<<side<<
1024       // pid type
1025       "isMuon="<<isMuon<<
1026       "isProton="<<isProton<<
1027       "isKaon="<<isKaon<<
1028       "isPion="<<isPion<<
1029       "isElectron="<<isElectron<<
1030       "isMIP="<<isMIP<<
1031       "isHighPt="<<isHighPt<<
1032       "mass="<<mass<<
1033       "dedxDef="<<dedxDef<<
1034       "dedxDefPion="<<dedxDefPion<<
1035       "dedxDefKaon="<<dedxDefKaon<<
1036       "dedxDefProton="<<dedxDefProton<<
1037       //
1038       "nclITS="<<nclITS<<
1039       "ncl="<<ncl<<
1040       "ncl21="<<ncl21<<
1041       "ncl20="<<ncl20<<
1042       //
1043       "ncl20All="<<ncl20All<<
1044       "ncl20OROC="<<ncl20OROC<<
1045       "ncl20IROC="<<ncl20IROC<<
1046       "ncl20OROC0="<<ncl20OROC0<<
1047       "ncl20OROC1="<<ncl20OROC1<<
1048       //
1049       "ncl20All4="<<ncl20All4<<
1050       "ncl20OROC4="<<ncl20OROC4<<
1051       "ncl20IROC4="<<ncl20IROC4<<
1052       "ncl20OROC04="<<ncl20OROC04<<
1053       "ncl20OROC14="<<ncl20OROC14<<
1054       //
1055       "ncl20All3="<<ncl20All3<<
1056       "ncl20OROC3="<<ncl20OROC3<<
1057       "ncl20IROC3="<<ncl20IROC3<<
1058       "ncl20OROC03="<<ncl20OROC03<<
1059       "ncl20OROC13="<<ncl20OROC13<<
1060       //
1061       "ncl21All="<<ncl21All<<
1062       "ncl21OROC="<<ncl21OROC<<
1063       "ncl21IROC="<<ncl21IROC<<
1064       "ncl21OROC0="<<ncl21OROC0<<
1065       "ncl21OROC1="<<ncl21OROC1<<
1066       //track properties
1067       "langle0="<<langle[0]<<
1068       "langle1="<<langle[1]<<
1069       "langle2="<<langle[2]<<
1070       "gangle0="<<gangle[0]<<   //global angle phi IROC
1071       "gangle1="<<gangle[1]<<   //                 OROC medium
1072       "gangle2="<<gangle[2]<<   //                 OROC long
1073       "L0="<<length[0]<<
1074       "L1="<<length[1]<<
1075       "L2="<<length[2]<<
1076       "p="<<p<<
1077       "pin="<<pin<<
1078       "pout="<<pout<<
1079       "tgl="<<tgl<<
1080       "track.="<<track<<
1081       "trackIn.="<<trackIn<<
1082       "trackOut.="<<trackOut<<
1083       "tpcOut.="<<tpcOut<<
1084       "seed.="<<seed<<
1085       "medianMIP0="<<medianMIP0<<    // median MIP position as estimated from the array of (not only) "MIPS"
1086       //dedx
1087       "truncUp.="<<&truncUp<<
1088       "truncDown.="<<&truncDown<<
1089       "vecAllMax.="<<&vecAllMax<<
1090       "vecIROCMax.="<<&vecIROCMax<<
1091       "vecOROCMax.="<<&vecOROCMax<<
1092       "vecOROC0Max.="<<&vecOROC0Max<<
1093       "vecOROC1Max.="<<&vecOROC1Max<<
1094       //
1095       "vecAllTot.="<<&vecAllTot<<
1096       "vecIROCTot.="<<&vecIROCTot<<
1097       "vecOROCTot.="<<&vecOROCTot<<
1098       "vecOROC0Tot.="<<&vecOROC0Tot<<
1099       "vecOROC1Tot.="<<&vecOROC1Tot<<
1100       //
1101       "vecAllTotLog.="<<&vecAllTotLog<<
1102       "vecIROCTotLog.="<<&vecIROCTotLog<<
1103       "vecOROCTotLog.="<<&vecOROCTotLog<<
1104       "vecOROC0TotLog.="<<&vecOROC0TotLog<<
1105       "vecOROC1TotLog.="<<&vecOROC1TotLog<<
1106       //
1107       "vecAllTotUp.="<<&vecAllTotUp<<
1108       "vecIROCTotUp.="<<&vecIROCTotUp<<
1109       "vecOROCTotUp.="<<&vecOROCTotUp<<
1110       "vecOROC0TotUp.="<<&vecOROC0TotUp<<
1111       "vecOROC1TotUp.="<<&vecOROC1TotUp<<
1112       //
1113       "vecAllTotDown.="<<&vecAllTotDown<<
1114       "vecIROCTotDown.="<<&vecIROCTotDown<<
1115       "vecOROCTotDown.="<<&vecOROCTotDown<<
1116       "vecOROC0TotDown.="<<&vecOROC0TotDown<<
1117       "vecOROC1TotDown.="<<&vecOROC1TotDown<<
1118       //
1119       "vecAllTotRMS.="<<&vecAllTotRMS<<
1120       "vecIROCTotRMS.="<<&vecIROCTotRMS<<
1121       "vecOROCTotRMS.="<<&vecOROCTotRMS<<
1122       "vecOROC0TotRMS.="<<&vecOROC0TotRMS<<
1123       "vecOROC1TotRMS.="<<&vecOROC1TotRMS<<
1124       //
1125       "vecAllTotM2.="<<&vecAllTotM2<<
1126       "vecIROCTotM2.="<<&vecIROCTotM2<<
1127       "vecOROCTotM2.="<<&vecOROCTotM2<<
1128       "vecOROC0TotM2.="<<&vecOROC0TotM2<<
1129       "vecOROC1TotM2.="<<&vecOROC1TotM2<<
1130       //
1131       "vecAllTotMS.="<<&vecAllTotMS<<
1132       "vecIROCTotMS.="<<&vecIROCTotMS<<
1133       "vecOROCTotMS.="<<&vecOROCTotMS<<
1134       "vecOROC0TotMS.="<<&vecOROC0TotMS<<
1135       "vecOROC1TotMS.="<<&vecOROC1TotMS<<
1136       "\n";
1137   }
1138 }
1139
1140
1141
1142
1143 void AliTPCcalibGainMult::ProcessV0s(AliESDEvent * event){
1144   //
1145   // Select the K0s and gamma  - and sign daughter products
1146   //
1147   TTreeSRedirector * pcstream =  GetDebugStreamer();
1148   AliKFParticle::SetField(event->GetMagneticField());
1149   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1150   if (!esdFriend) {
1151     //Printf("ERROR: esdFriend not available");
1152    return;
1153   }
1154   if (esdFriend->TestSkipBit()) return;
1155   //
1156   //
1157   static const TDatabasePDG *pdg = TDatabasePDG::Instance();
1158   const Double_t kChi2Cut=5;
1159   const Double_t kMinR=2;
1160   const Int_t    kMinNcl=110;
1161   const Double_t kMinREl=5;
1162   const Double_t kMaxREl=70;
1163   //
1164   Int_t nv0 = event->GetNumberOfV0s();
1165   AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex();
1166   AliKFVertex kfvertex=*vertex;
1167   //
1168   for (Int_t iv0=0;iv0<nv0;iv0++){
1169     AliESDv0 *v0 = event->GetV0(iv0);
1170     if (!v0) continue;
1171     if (v0->GetOnFlyStatus()<0.5) continue;
1172     if (v0->GetPindex()<0) continue;
1173     if (v0->GetNindex()<0) continue;
1174     if (TMath::Max(v0->GetPindex(), v0->GetNindex())>event->GetNumberOfTracks()) continue;
1175     //
1176     //
1177     AliExternalTrackParam pp=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamP())):(*(v0->GetParamN()));
1178     AliExternalTrackParam pn=(v0->GetParamP()->GetSign()>0) ? (*(v0->GetParamN())):(*(v0->GetParamP()));
1179     AliKFParticle kfp1( pp, 211 );
1180     AliKFParticle kfp2( pn, -211 );
1181     //
1182     AliKFParticle *v0KFK0 = new AliKFParticle(kfp1,kfp2);
1183     AliKFParticle *v0KFK0CV = new AliKFParticle(*v0KFK0);
1184     v0KFK0CV->SetProductionVertex(kfvertex);
1185     v0KFK0CV->TransportToProductionVertex();
1186     Double_t chi2K0 = v0KFK0CV->GetChi2();
1187     if (chi2K0>kChi2Cut) continue;
1188     if (v0->GetRr()<kMinR) continue;
1189     Bool_t isOKC=TMath::Max(v0->GetCausalityP()[0],v0->GetCausalityP()[1])<0.7&&TMath::Min(v0->GetCausalityP()[2],v0->GetCausalityP()[3])>0.2;
1190     //
1191     Double_t effMass22=v0->GetEffMass(2,2);
1192     Double_t effMass42=v0->GetEffMass(4,2);
1193     Double_t effMass24=v0->GetEffMass(2,4);
1194     Double_t effMass00=v0->GetEffMass(0,0);
1195     AliKFParticle *v0KFK0CVM = new AliKFParticle(*v0KFK0CV);
1196     v0KFK0CVM->SetMassConstraint(pdg->GetParticle("K_S0")->Mass());
1197     Bool_t isV0= kFALSE;
1198     //
1199     Double_t d22   = TMath::Abs(effMass22-pdg->GetParticle("K_S0")->Mass());
1200     Double_t d42   = TMath::Abs(effMass42-pdg->GetParticle("Lambda0")->Mass());
1201     Double_t d24   = TMath::Abs(effMass24-pdg->GetParticle("Lambda0")->Mass());
1202     Double_t d00   = TMath::Abs(effMass00);
1203     //
1204     Bool_t isKaon      = d22<0.01 && d22< 0.3 * TMath::Min(TMath::Min(d42,d24),d00);
1205     Bool_t isLambda    = d42<0.01 && d42< 0.3 * TMath::Min(TMath::Min(d22,d24),d00);
1206     Bool_t isAntiLambda= d24<0.01 && d24< 0.3 * TMath::Min(TMath::Min(d22,d42),d00);
1207     Bool_t isGamma     = d00<0.02 && d00< 0.3 * TMath::Min(TMath::Min(d42,d24),d22);
1208     //
1209     if (isGamma  &&  (isKaon||isLambda||isAntiLambda)) continue;
1210     if (isLambda &&  (isKaon||isGamma||isAntiLambda)) continue;
1211     if (isKaon   &&  (isLambda||isGamma||isAntiLambda)) continue;
1212     Double_t sign= v0->GetParamP()->GetSign()* v0->GetParamN()->GetSign();
1213     if (sign>0) continue;
1214     isV0=isKaon||isLambda||isAntiLambda||isGamma;
1215     if (!(isV0)) continue;
1216     if (isGamma&&v0->GetRr()<kMinREl) continue;
1217     if (isGamma&&v0->GetRr()>kMaxREl) continue;
1218     if (!isOKC) continue;
1219     //
1220     Int_t pindex = (v0->GetParamP()->GetSign()>0) ? v0->GetPindex() : v0->GetNindex();
1221     Int_t nindex = (v0->GetParamP()->GetSign()>0) ? v0->GetNindex() : v0->GetPindex();
1222     AliESDtrack * trackP = event->GetTrack(pindex);
1223     AliESDtrack * trackN = event->GetTrack(nindex);
1224     if (!trackN) continue;
1225     if (!trackP) continue;
1226     Int_t nclP= (Int_t)trackP->GetTPCClusterInfo(2,1);
1227     Int_t nclN= (Int_t)trackN->GetTPCClusterInfo(2,1);
1228     if (TMath::Min(nclP,nclN)<kMinNcl) continue;
1229     Double_t eta = TMath::Max(TMath::Abs(trackP->Eta()), TMath::Abs(trackN->Eta()));
1230     if (TMath::Abs(eta)>1) continue;
1231     //
1232     //
1233     AliESDfriendTrack *friendTrackP = esdFriend->GetTrack(pindex);
1234     AliESDfriendTrack *friendTrackN = esdFriend->GetTrack(nindex);
1235     if (!friendTrackP) continue;
1236     if (!friendTrackN) continue;
1237     TObject *calibObject;
1238     AliTPCseed *seedP = 0;
1239     AliTPCseed *seedN = 0;
1240     for (Int_t l=0;(calibObject=friendTrackP->GetCalibObject(l));++l) {
1241       if ((seedP=dynamic_cast<AliTPCseed*>(calibObject))) break;
1242     }
1243     for (Int_t l=0;(calibObject=friendTrackN->GetCalibObject(l));++l) {
1244       if ((seedN=dynamic_cast<AliTPCseed*>(calibObject))) break;
1245     }
1246     if (isGamma){
1247       if ( TMath::Abs((trackP->GetTPCsignal()/(trackN->GetTPCsignal()+0.0001)-1)>0.3)) continue;
1248     }
1249     if (isGamma)   (*fPIDMatrix)(pindex, 0)+=2;
1250     if (isGamma)   (*fPIDMatrix)(nindex, 0)+=2;
1251     //
1252     if (isKaon)    (*fPIDMatrix)(pindex, 2)+=2;
1253     if (isKaon)    (*fPIDMatrix)(nindex, 2)+=2;
1254     //
1255     //
1256     if (pcstream){
1257       (*pcstream)<<"v0s"<<
1258         "isGamma="<<isGamma<<
1259         "isKaon="<<isKaon<<
1260         "isLambda="<<isLambda<<
1261         "isAntiLambda="<<isAntiLambda<<
1262         "chi2="<<chi2K0<<
1263         "trackP.="<<trackP<<
1264         "trackN.="<<trackN<<
1265         "v0.="<<v0<<
1266         "\n";
1267     }
1268   }
1269 }
1270
1271
1272
1273
1274 void AliTPCcalibGainMult::ProcessCosmic(const AliESDEvent * event) {
1275   //
1276   // Find cosmic pairs trigger by random trigger
1277   //
1278   //
1279   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1280   AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
1281
1282   AliESDVertex *vertexSPD =  (AliESDVertex *)event->GetPrimaryVertexSPD();
1283   AliESDVertex *vertexTPC =  (AliESDVertex *)event->GetPrimaryVertexTPC();
1284   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1285   const Double_t kMinPt=4;
1286   const Double_t kMinPtMax=0.8;
1287   const Double_t kMinNcl=159*0.5;
1288   const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
1289   Int_t ntracks=event->GetNumberOfTracks();
1290   const Double_t kMaxImpact=80;
1291   //  Float_t dcaTPC[2]={0,0};
1292   // Float_t covTPC[3]={0,0,0};
1293
1294   UInt_t specie = event->GetEventSpecie();  // skip laser events
1295   if (specie==AliRecoParam::kCalib) return;
1296
1297
1298   for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
1299     AliESDtrack *track0 = event->GetTrack(itrack0);
1300     if (!track0) continue;
1301     if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
1302
1303     if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()<kMinPt) continue;
1304     if (track0->GetTPCncls()<kMinNcl) continue;
1305     if (TMath::Abs(track0->GetY())<2*kMaxDelta[0]) continue;
1306     if (TMath::Abs(track0->GetY())>kMaxImpact) continue;
1307     if (track0->GetKinkIndex(0)>0) continue;
1308     const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
1309     //rm primaries
1310     //
1311     for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
1312       AliESDtrack *track1 = event->GetTrack(itrack1);
1313       if (!track1) continue;
1314       if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
1315       if (track1->GetKinkIndex(0)>0) continue;
1316       if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()<kMinPt) continue;
1317       if (track1->GetTPCncls()<kMinNcl) continue;
1318       if (TMath::Abs(AliTracker::GetBz())>1&&TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
1319       if (TMath::Abs(track1->GetY())<2*kMaxDelta[0]) continue;
1320       if (TMath::Abs(track1->GetY())>kMaxImpact) continue;
1321       //
1322       const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
1323       //
1324       Bool_t isPair=kTRUE;
1325       for (Int_t ipar=0; ipar<5; ipar++){
1326         if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
1327         if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
1328       }
1329       if (!isPair) continue;
1330       if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
1331       //delta with correct sign
1332       if  (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y   opposite sign
1333       if  (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
1334       if  (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
1335       if (!isPair) continue;
1336       TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1337       Int_t eventNumber = event->GetEventNumberInFile();
1338       Bool_t hasFriend=(esdFriend) ? (esdFriend->GetTrack(itrack0)!=0):0;
1339       Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);
1340       printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS);
1341       //
1342       //
1343       TTreeSRedirector * pcstream =  GetDebugStreamer();
1344       Int_t ntracksSPD = vertexSPD->GetNContributors();
1345       Int_t ntracksTPC = vertexTPC->GetNContributors();
1346       //
1347       if (pcstream){
1348         (*pcstream)<<"cosmicPairsAll"<<
1349           "run="<<fRun<<              //  run number
1350           "event="<<fEvent<<          //  event number
1351           "time="<<fTime<<            //  time stamp of event
1352           "trigger="<<fTrigger<<      //  trigger
1353           "triggerClass="<<&fTriggerClass<<      //  trigger
1354           "bz="<<fMagF<<             //  magnetic field
1355           //
1356           "nSPD="<<ntracksSPD<<
1357           "nTPC="<<ntracksTPC<<
1358           "vSPD.="<<vertexSPD<<         //primary vertex -SPD
1359           "vTPC.="<<vertexTPC<<         //primary vertex -TPC
1360           "t0.="<<track0<<              //track0
1361           "t1.="<<track1<<              //track1
1362           "\n";
1363       }
1364       //
1365       AliESDfriendTrack *friendTrack0 = esdFriend->GetTrack(itrack0);
1366       if (!friendTrack0) continue;
1367       AliESDfriendTrack *friendTrack1 = esdFriend->GetTrack(itrack1);
1368       if (!friendTrack1) continue;
1369       TObject *calibObject;
1370       AliTPCseed *seed0 = 0;
1371       AliTPCseed *seed1 = 0;
1372       //
1373       for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) {
1374         if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1375       }
1376       for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) {
1377         if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1378       }
1379       //
1380       if (pcstream){
1381         (*pcstream)<<"cosmicPairs"<<
1382           "run="<<fRun<<              //  run number
1383           "event="<<fEvent<<          //  event number
1384           "time="<<fTime<<            //  time stamp of event
1385           "trigger="<<fTrigger<<      //  trigger
1386           "triggerClass="<<&fTriggerClass<<      //  trigger
1387           "bz="<<fMagF<<             //  magnetic field
1388           //
1389           "nSPD="<<ntracksSPD<<
1390           "nTPC="<<ntracksTPC<<
1391           "vSPD.="<<vertexSPD<<         //primary vertex -SPD
1392           "vTPC.="<<vertexTPC<<         //primary vertex -TPC
1393           "t0.="<<track0<<              //track0
1394           "t1.="<<track1<<              //track1
1395           "ft0.="<<friendTrack0<<       //track0
1396           "ft1.="<<friendTrack1<<       //track1
1397           "s0.="<<seed0<<               //track0
1398           "s1.="<<seed1<<               //track1
1399           "\n";
1400       }
1401       if (!seed0) continue;
1402       if (!seed1) continue;
1403       Int_t nclA0=0, nclC0=0;     // number of clusters
1404       Int_t nclA1=0, nclC1=0;     // number of clusters
1405       //
1406       for (Int_t irow=0; irow<159; irow++){
1407         AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1408         AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1409         if (cluster0){
1410           if (cluster0->GetQ()>0 &&  cluster0->GetDetector()%36<18)  nclA0++;
1411           if (cluster0->GetQ()>0 &&  cluster0->GetDetector()%36>=18) nclC0++;
1412         }
1413         if (cluster1){
1414           if (cluster1->GetQ()>0 &&  cluster1->GetDetector()%36<18)  nclA1++;
1415           if (cluster1->GetQ()>0 &&  cluster1->GetDetector()%36>=18) nclC1++;
1416         }
1417       }
1418       Int_t cosmicType=0;  // types of cosmic topology
1419       if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side
1420       if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side
1421       if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side
1422       if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side
1423       if (cosmicType<2) continue; // use only crossing tracks
1424       //
1425       Double_t deltaTimeCluster=0;
1426       deltaTimeCluster=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth();
1427       if (nclA0>nclC0) deltaTimeCluster*=-1; // if A side track
1428       //
1429       for (Int_t irow=0; irow<159; irow++){
1430         AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1431         if (cluster0 &&cluster0->GetX()>10){
1432           Double_t x0[3]={cluster0->GetRow(),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster};
1433           Int_t index0[1]={cluster0->GetDetector()};
1434           transform->Transform(x0,index0,0,1);
1435           cluster0->SetX(x0[0]);
1436           cluster0->SetY(x0[1]);
1437           cluster0->SetZ(x0[2]);
1438           //
1439         }
1440         AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1441         if (cluster1&&cluster1->GetX()>10){
1442           Double_t x1[3]={cluster1->GetRow(),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster};
1443           Int_t index1[1]={cluster1->GetDetector()};
1444           transform->Transform(x1,index1,0,1);
1445           cluster1->SetX(x1[0]);
1446           cluster1->SetY(x1[1]);
1447           cluster1->SetZ(x1[2]);
1448         }
1449       }
1450       //
1451       //
1452       if (fPIDMatrix){
1453         (*fPIDMatrix)(itrack0,1)+=4;  //
1454         (*fPIDMatrix)(itrack1,1)+=4;  //
1455       }
1456     }
1457   }
1458 }
1459
1460
1461
1462 void AliTPCcalibGainMult::ProcessKinks(const AliESDEvent * event){
1463   //
1464   //
1465   //
1466   AliKFParticle::SetField(event->GetMagneticField());
1467   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1468   if (!esdFriend) {
1469     //Printf("ERROR: esdFriend not available");
1470     return;
1471   }
1472   //  if (esdFriend->TestSkipBit()) return;
1473   //
1474   //
1475   const Double_t kChi2Cut=10;
1476   const Double_t kMinR=100;
1477   const Double_t kMaxR=230;
1478   const Int_t    kMinNcl=110;
1479   //
1480   Int_t nkinks = event->GetNumberOfKinks();
1481   AliESDVertex *vertex= (AliESDVertex *)event->GetPrimaryVertex();
1482   AliKFVertex kfvertex=*vertex;
1483   TTreeSRedirector * pcstream =  GetDebugStreamer();
1484   //
1485   for (Int_t ikink=0;ikink<nkinks;ikink++){
1486     AliESDkink *kink = event->GetKink(ikink);
1487     if (!kink) continue;
1488     if (kink->GetIndex(0)<0) continue;
1489     if (kink->GetIndex(1)<0) continue;
1490     if (TMath::Max(kink->GetIndex(1), kink->GetIndex(0))>event->GetNumberOfTracks()) continue;
1491     //
1492     //
1493     AliExternalTrackParam pd=kink->RefParamDaughter();
1494     AliExternalTrackParam pm=kink->RefParamMother();
1495     AliKFParticle kfpd( pd, 211 );
1496     AliKFParticle kfpm( pm, -13 );
1497     //
1498     AliKFParticle *v0KF = new AliKFParticle(kfpm,kfpd);
1499     v0KF->SetVtxGuess(kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]);
1500     Double_t chi2 = v0KF->GetChi2();
1501     AliESDtrack * trackM = event->GetTrack(kink->GetIndex(0));
1502     AliESDtrack * trackD = event->GetTrack(kink->GetIndex(1));
1503     if (!trackM) continue;
1504     if (!trackD) continue;
1505     Int_t nclM= (Int_t)trackM->GetTPCClusterInfo(2,1);
1506     Int_t nclD= (Int_t)trackD->GetTPCClusterInfo(2,1);
1507     Double_t eta = TMath::Max(TMath::Abs(trackM->Eta()), TMath::Abs(trackD->Eta()));
1508     Double_t kx= v0KF->GetX();
1509     Double_t ky= v0KF->GetY();
1510     Double_t kz= v0KF->GetZ();
1511     Double_t ex= v0KF->GetErrX();
1512     Double_t ey= v0KF->GetErrY();
1513     Double_t ez= v0KF->GetErrZ();
1514     //
1515     Double_t radius=TMath::Sqrt(kx*kx+ky*ky);
1516     Double_t alpha=TMath::ATan2(ky,kx);
1517     if (!pd.Rotate(alpha)) continue;
1518     if (!pm.Rotate(alpha)) continue;
1519     if (!pd.PropagateTo(radius,event->GetMagneticField())) continue;
1520     if (!pm.PropagateTo(radius,event->GetMagneticField())) continue;
1521     Double_t pos[2]={0,kz};
1522     Double_t cov[3]={ex*ex+ey*ey,0,ez*ez};
1523     pd.Update(pos,cov);
1524     pm.Update(pos,cov);
1525     //
1526     if (pcstream){
1527       (*pcstream)<<"kinks"<<
1528         "eta="<<eta<<
1529         "nclM="<<nclM<<
1530         "nclD="<<nclD<<
1531         "kink.="<<kink<<
1532         "trackM.="<<trackM<<
1533         "trackD.="<<trackD<<
1534         "pm.="<<&pm<<             //updated parameters
1535         "pd.="<<&pd<<             // updated parameters
1536         "v0KF.="<<v0KF<<
1537         "chi2="<<chi2<<
1538         "\n";
1539     }
1540     /*
1541       TCut cutQ="chi2<10&&kink.fRr>90&&kink.fRr<220";
1542       TCut cutRD="20*sqrt(pd.fC[14])<abs(pd.fP[4])&&trackD.fTPCsignal>10&&trackD.fTPCsignalN>50";
1543
1544     */
1545     //
1546     if (chi2>kChi2Cut) continue;
1547     if (kink->GetR()<kMinR) continue;
1548     if (kink->GetR()>kMaxR) continue;
1549     if ((nclM+nclD)<kMinNcl) continue;
1550     if (TMath::Abs(eta)>1) continue;
1551     //
1552     //
1553     AliESDfriendTrack *friendTrackM = esdFriend->GetTrack(kink->GetIndex(0));
1554     AliESDfriendTrack *friendTrackD = esdFriend->GetTrack(kink->GetIndex(1));
1555     if (!friendTrackM) continue;
1556     if (!friendTrackD) continue;
1557     TObject *calibObject;
1558     AliTPCseed *seedM = 0;
1559     AliTPCseed *seedD = 0;
1560     for (Int_t l=0;(calibObject=friendTrackM->GetCalibObject(l));++l) {
1561       if ((seedM=dynamic_cast<AliTPCseed*>(calibObject))) break;
1562     }
1563     for (Int_t l=0;(calibObject=friendTrackD->GetCalibObject(l));++l) {
1564       if ((seedD=dynamic_cast<AliTPCseed*>(calibObject))) break;
1565     }
1566   }
1567 }
1568
1569 void AliTPCcalibGainMult::DumpHPT(const AliESDEvent * event){
1570   //
1571   // Function to select the HPT tracks and events
1572   // It is used to select event with HPT - list used later for the raw data downloading
1573   //                                     - and reconstruction
1574   // Not actualy used for the calibration of the data
1575
1576   TTreeSRedirector * pcstream =  GetDebugStreamer();
1577   AliKFParticle::SetField(event->GetMagneticField());
1578   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1579   if (!esdFriend) {
1580     //Printf("ERROR: esdFriend not available");
1581    return;
1582   }
1583   if (esdFriend->TestSkipBit()) return;
1584
1585   Int_t ntracks=event->GetNumberOfTracks();
1586   //
1587   for (Int_t i=0;i<ntracks;++i) {
1588     //
1589     AliESDtrack *track = event->GetTrack(i);
1590     if (!track) continue;
1591     if (track->Pt()<4) continue;
1592     UInt_t status = track->GetStatus();
1593     //
1594     AliExternalTrackParam * trackIn  = (AliExternalTrackParam *)track->GetInnerParam();
1595     if (!trackIn) continue;
1596     if ((status&AliESDtrack::kTPCrefit)==0) continue;
1597     if ((status&AliESDtrack::kITSrefit)==0) continue;
1598     AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
1599     if (!friendTrack) continue;
1600     AliExternalTrackParam * itsOut = (AliExternalTrackParam *)(friendTrack->GetITSOut());
1601     if (!itsOut) continue;
1602     AliExternalTrackParam * itsOut2 = (AliExternalTrackParam *)(friendTrack->GetITSOut()->Clone());
1603     AliExternalTrackParam * tpcIn2 = (AliExternalTrackParam *)(trackIn->Clone());
1604     if (!itsOut2->Rotate(trackIn->GetAlpha())) continue;
1605     //Double_t xmiddle=0.5*(itsOut2->GetX()+tpcIn2->GetX());
1606     Double_t xmiddle=(itsOut2->GetX());
1607     if (!itsOut2->PropagateTo(xmiddle,event->GetMagneticField())) continue;
1608     if (!tpcIn2->PropagateTo(xmiddle,event->GetMagneticField())) continue;
1609     //
1610     AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
1611     if (!tpcInner) continue;
1612     tpcInner->Rotate(track->GetAlpha());
1613     tpcInner->PropagateTo(track->GetX(),event->GetMagneticField());
1614     //
1615     // tpc constrained
1616     //
1617     AliExternalTrackParam * tpcInnerC = (AliExternalTrackParam *)(track->GetTPCInnerParam()->Clone());
1618     if (!tpcInnerC) continue;
1619     tpcInnerC->Rotate(track->GetAlpha());
1620     tpcInnerC->PropagateTo(track->GetX(),event->GetMagneticField());
1621     Double_t dz[2],cov[3];
1622     AliESDVertex *vtx= (AliESDVertex *)event->GetPrimaryVertex();
1623
1624     if (!tpcInnerC->PropagateToDCA(vtx, event->GetMagneticField(), 3, dz, cov)) continue;
1625     Double_t covar[6]; vtx->GetCovMatrix(covar);
1626     Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
1627     Double_t c[3]={covar[2],0.,covar[5]};
1628     //
1629     Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
1630     tpcInnerC->Update(p,c);
1631
1632     if (pcstream){
1633       (*pcstream)<<"hpt"<<
1634         "run="<<fRun<<
1635         "time="<<fTime<<
1636         "vertex="<<vtx<<
1637         "bz="<<fMagF<<
1638         "track.="<<track<<
1639         "tpcInner.="<<tpcInner<<
1640         "tpcInnerC.="<<tpcInnerC<<
1641         "chi2C="<<chi2C<<
1642         //
1643         "its.="<<itsOut<<
1644         "its2.="<<itsOut2<<
1645         "tpc.="<<trackIn<<
1646         "tpc2.="<<tpcIn2<<
1647         "\n";
1648     }
1649   }
1650 }
1651
1652
1653
1654 void AliTPCcalibGainMult::ProcessTOF(const AliESDEvent * event){
1655   //
1656   // 1. Loop over tracks
1657   // 2. Fit T0
1658   // 3. Sign positivelly identified tracks
1659   //
1660   const Double_t kMaxDelta=1000;
1661   const Double_t kOrbit=50000; // distance in the time beween two orbits in the TOF time units  - 50000=50 ns
1662   const Double_t kMaxD=20;
1663   const Double_t kRMS0=200;
1664   const Double_t kMaxDCAZ=10;
1665   AliESDVertex *vtx= (AliESDVertex *)event->GetPrimaryVertex();
1666   //
1667   TTreeSRedirector * pcstream =  GetDebugStreamer();
1668   AliKFParticle::SetField(event->GetMagneticField());
1669   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
1670   if (!esdFriend) {
1671     //Printf("ERROR: esdFriend not available");
1672    return;
1673   }
1674   //if (esdFriend->TestSkipBit()) return;
1675
1676   Int_t ntracks=event->GetNumberOfTracks();
1677   //
1678   Double_t deltaTPion[10000];
1679   Double_t medianT0=0;
1680   Double_t meanT0=0;
1681   Double_t rms=10000;
1682   Int_t counter=0;
1683   //
1684   // Get Median time for pion hypothesy
1685   //
1686   for (Int_t iter=0; iter<3; iter++){
1687     counter=0;
1688     for (Int_t i=0;i<ntracks;++i) {
1689       //
1690       AliESDtrack *track = event->GetTrack(i);
1691       if (!track) continue;
1692       if (!track->IsOn(AliESDtrack::kTIME)) continue;
1693       if (TMath::Abs(track->GetZ())>kMaxDCAZ) continue;         // remove overlaped events
1694       if (TMath::Abs(track->GetTOFsignalDz())>kMaxD) continue;
1695       Double_t times[1000];
1696       track->GetIntegratedTimes(times);
1697       Int_t norbit=TMath::Nint((track->GetTOFsignal()-times[2])/kOrbit);
1698       Double_t torbit=norbit*kOrbit;
1699       if (iter==1 &&TMath::Abs(times[2]-times[3])<3*rms) continue;  // skip umbigous points - kaon pion
1700       //
1701       Int_t indexBest=2;
1702       if (iter>1){
1703         for (Int_t j=3; j<5; j++)
1704           if (TMath::Abs(track->GetTOFsignal()-times[j]-torbit-medianT0)<TMath::Abs(track->GetTOFsignal()-times[j]-torbit-medianT0)) indexBest=j;
1705       }
1706       //
1707       if (iter>0) if (TMath::Abs(track->GetTOFsignal()-times[indexBest]-torbit-medianT0)>3*(kRMS0+rms)) continue;
1708       if (iter>0) if (TMath::Abs(track->GetTOFsignal()-times[indexBest]-torbit-medianT0)>kMaxDelta) continue;
1709       deltaTPion[counter]=track->GetTOFsignal()-times[indexBest]-torbit;
1710       counter++;
1711     }
1712     if (counter<2) return;
1713     medianT0=TMath::Median(counter,deltaTPion);
1714     meanT0=TMath::Median(counter,deltaTPion);
1715     rms=TMath::RMS(counter,deltaTPion);
1716   }
1717   if (counter<3) return;
1718   //
1719   // Dump
1720   //
1721   for (Int_t i=0;i<ntracks;++i) {
1722     //
1723     AliESDtrack *track = event->GetTrack(i);
1724     if (!track) continue;
1725     if (!track->IsOn(AliESDtrack::kTIME)) continue;
1726     if (TMath::Abs(track->GetZ())>kMaxDCAZ) continue;          //remove overlapped events
1727     if (TMath::Abs(track->GetTOFsignalDz())>kMaxD) continue;
1728     Double_t times[1000];
1729     track->GetIntegratedTimes(times);
1730     Int_t norbit=TMath::Nint((track->GetTOFsignal()-times[2])/kOrbit);
1731     Double_t torbit=norbit*kOrbit;
1732     if (rms<=0) continue;
1733     //
1734     Double_t tPion  = (track->GetTOFsignal()-times[2]-medianT0-torbit);
1735     Double_t tKaon  = (track->GetTOFsignal()-times[3]-medianT0-torbit);
1736     Double_t tProton= (track->GetTOFsignal()-times[4]-medianT0-torbit);
1737     Double_t tElectron= (track->GetTOFsignal()-times[0]-medianT0-torbit);
1738     //
1739     Bool_t isPion   = (TMath::Abs(tPion/rms)<6)   && TMath::Abs(tPion)<(TMath::Min(TMath::Abs(tKaon), TMath::Abs(tProton))-rms);
1740     Bool_t isKaon   = (TMath::Abs(tKaon/rms)<3)   && TMath::Abs(tKaon)<0.2*(TMath::Min(TMath::Abs(tPion), TMath::Abs(tProton))-3*rms);
1741     Bool_t isProton = (TMath::Abs(tProton/rms)<6) && TMath::Abs(tProton)<0.5*(TMath::Min(TMath::Abs(tKaon), TMath::Abs(tPion))-rms);
1742     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);
1743
1744     if (isPion)   (*fPIDMatrix)(i,2)+=1;
1745     if (isKaon)   (*fPIDMatrix)(i,3)+=1;
1746     if (isProton) (*fPIDMatrix)(i,4)+=1;
1747     //    if (isElectron) (*fPIDMatrix)(i,0)+=1;
1748     //
1749     if (pcstream){
1750       // debug streamer to dump the information
1751     (*pcstream)<<"tof"<<
1752       "isPion="<<isPion<<
1753       "isKaon="<<isKaon<<
1754       "isProton="<<isProton<<
1755       "isElectron="<<isElectron<<
1756       //
1757       "counter="<<counter<<
1758       "torbit="<<torbit<<
1759       "norbit="<<norbit<<
1760       "medianT0="<<medianT0<<
1761       "meanT0="<<meanT0<<
1762       "rmsT0="<<rms<<
1763       "track.="<<track<<
1764       "vtx.="<<vtx<<
1765       "\n";
1766     }
1767
1768   }
1769   /*
1770     tof->SetAlias("isProton","(abs(track.fTOFsignal-track.fTrackTime[4]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)-rmsT0))");
1771     tof->SetAlias("isPion","(abs(track.fTOFsignal-track.fTrackTime[2]-medianT0-torbit)<(0.5*abs(track.fTOFsignal-track.fTrackTime[3]-medianT0-torbit)-rmsT0))");
1772     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))");
1773
1774    */
1775
1776 }
1777
1778
1779 // void AliTPCcalibGainMult::Terminate(){
1780 //   //
1781 //   // Terminate function
1782 //   // call base terminate + Eval of fitters
1783 //   //
1784 //    Info("AliTPCcalibGainMult","Terminate");
1785 //    TTreeSRedirector *pcstream = GetDebugStreamer();
1786 //    if (pcstream){
1787 //      TTreeStream &stream = (*pcstream)<<"dump";
1788 //      TTree* tree = stream.GetTree();
1789 //      if (tree) if ( tree->GetEntries()>0){
1790 //        TObjArray *array =  tree->GetListOfBranches();
1791 //        for (Int_t i=0; i<array->GetEntries(); i++) {TBranch * br = (TBranch *)array->At(i); br->SetAddress(0);}
1792 //        gDirectory=gROOT;
1793 //        fdEdxTree=tree->CloneTree(10000);
1794 //      }
1795 //    }
1796 //    AliTPCcalibBase::Terminate();
1797 // }
1798