TPC/Calib/AliTPCcalibGainMult.cxx

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