TPC/AliTPCcalibGainMult.cxx

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