TPC/AliTPCcalibGainMult.cxx

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