TPC/Calib/AliTPCcalibGainMult.cxx

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