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