TPC/AliTPCcalibGainMult.cxx

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