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