PWG1/AliMCTrackingTestTask.cxx

   1 //
   2 // This class is the task to check the
   3 // Propagation and Update method used in the
   4 //               1. AliExternalTrackParam
   5 //               2. AliTracker
   6 //
   7 // Pure Monte-Carlo data used, not influence of detectors
   8
   9 // Input - TParticle + Array of track references - (Points alogn track trajectories)
  10 // Output - Trees with track references - no histograms
  11 //          MC tree -  test for material budget correction
  12 //                     see function ProcessRefTracker
  13 //          MCupdate tree - test for correctness of propagation and update
  14 //                     see function AliMCTrackingTestTask::FitTrackRefs
  15 //
  16 // Principle - Creates AliExternalTrackParam form 1 Track Refernece -
  17 //             Propagate it to other
  18 // Magnetic field and the geometry has to be created before using it
  19
  20 //
  21 //
  22 //
  23
  24 //
  25 // ROOT includes
  26 #include <TChain.h>
  27 #include <TMath.h>
  28 #include <TVectorD.h>
  29 #include <TSystem.h>
  30 #include <TFile.h>
  31
  32 // ALIROOT includes
  33 #include <TTreeStream.h>
  34 #include <AliAnalysisManager.h>
  35 #include <AliESDInputHandler.h>
  36 #include "AliStack.h"
  37 #include "AliMCEvent.h"
  38 #include "AliMCEventHandler.h"
  39
  40 #include <AliESD.h>
  41 #include "AliMCTrackingTestTask.h"
  42 #include "AliGenInfoMaker.h"
  43 #include "AliHelix.h"
  44
  45 //
  46 #include "AliMCInfo.h"
  47 #include "AliComparisonObject.h"
  48 #include "AliESDRecInfo.h"
  49 #include "AliTPCParamSR.h"
  50 #include "AliTracker.h"
  51 #include "AliTPCseed.h"
  52
  53 // STL includes
  54 #include <iostream>
  55
  56 using namespace std;
  57
  58 ClassImp(AliMCTrackingTestTask)
  59
  60 //________________________________________________________________________
  61 AliMCTrackingTestTask::AliMCTrackingTestTask() :
  62   AliAnalysisTask(),
  63   fMCinfo(0),     //! MC event handler
  64   fESD(0),
  65   fDebugStreamer(0),
  66   fStreamLevel(0),
  67   fDebugLevel(0),
  68   fDebugOutputPath()
  69 {
  70   //
  71   // Default constructor (should not be used)
  72   //
  73 }
  74
  75 AliMCTrackingTestTask::AliMCTrackingTestTask(const AliMCTrackingTestTask& /*info*/) :
  76   AliAnalysisTask(),
  77   fMCinfo(0),     //! MC event handler
  78   fESD(0),
  79   //
  80   fDebugStreamer(0),
  81   fStreamLevel(0),
  82   fDebugLevel(),
  83   fDebugOutputPath()
  84 {
  85   //
  86   // Default constructor
  87   //
  88 }
  89
  90
  91
  92 //________________________________________________________________________
  93 AliMCTrackingTestTask::AliMCTrackingTestTask(const char *name) :
  94   AliAnalysisTask(name, "AliMCTrackingTestTask"),
  95   fMCinfo(0),     //! MC event handler
  96   fESD(0),
  97   fDebugStreamer(0),
  98   fStreamLevel(0),
  99   fDebugLevel(0),
 100   fDebugOutputPath()
 101 {
 102   //
 103   // Normal constructor
 104   //
 105   // Input slot #0 works with a TChain
 106   DefineInput(0, TChain::Class());
 107   // Output slot #0 writes into a TList
 108   DefineOutput(0, TObjArray::Class());
 109   //
 110   //
 111 }
 112
 113 AliMCTrackingTestTask::~AliMCTrackingTestTask(){
 114   //
 115   //
 116   //
 117   if (fDebugLevel>0)  printf("AliMCTrackingTestTask::~AliMCTrackingTestTask\n");
 118   if (fDebugStreamer) delete fDebugStreamer;
 119   fDebugStreamer=0;
 120 }
 121
 122
 123 //________________________________________________________________________
 124 void AliMCTrackingTestTask::ConnectInputData(Option_t *)
 125 {
 126   //
 127   // Connect the input data
 128   //
 129   if(fDebugLevel>3)
 130     cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
 131
 132   TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
 133   if (!tree) {
 134     //Printf("ERROR: Could not read chain from input slot 0");
 135   }
 136   else {
 137     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 138     if (!esdH) {
 139       //Printf("ERROR: Could not get ESDInputHandler");
 140     }
 141     else {
 142       fESD = esdH->GetEvent();
 143       //Printf("*** CONNECTED NEW EVENT ****");
 144     }
 145   }
 146   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 147   mcinfo->SetReadTR(kTRUE);
 148
 149   fMCinfo = mcinfo->MCEvent();
 150
 151
 152 }
 153
 154
 155
 156
 157
 158
 159 //________________________________________________________________________
 160 void AliMCTrackingTestTask::CreateOutputObjects()
 161 {
 162   //
 163   // Connect the output objects
 164   //
 165   if(fDebugLevel>3)
 166     cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
 167
 168 }
 169
 170
 171 //________________________________________________________________________
 172 void AliMCTrackingTestTask::Exec(Option_t *) {
 173   //
 174   // Execute analysis for current event
 175   //
 176
 177   if(fDebugLevel>3)
 178     cout << "AliMCTrackingTestTask::Exec()" << endl;
 179
 180
 181   // If MC has been connected
 182
 183   if (!fMCinfo){
 184     cout << "Not MC info\n" << endl;
 185   }else{
 186     ProcessMCInfo();
 187     //mcinfo->Print();
 188     //DumpInfo();
 189   }
 190   //
 191 }
 192
 193
 194
 195
 196 //________________________________________________________________________
 197 void AliMCTrackingTestTask::Terminate(Option_t *) {
 198     //
 199     // Terminate loop
 200     //
 201   if(fDebugLevel>3)
 202     printf("AliMCTrackingTestTask: Terminate() \n");
 203   //
 204   if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
 205   if (fDebugStreamer) delete fDebugStreamer;
 206   fDebugStreamer = 0;
 207   return;
 208 }
 209
 210
 211
 212 TTreeSRedirector *AliMCTrackingTestTask::GetDebugStreamer(){
 213   //
 214   // Get Debug streamer
 215   // In case debug streamer not yet initialized and StreamLevel>0 create new one
 216   //
 217   if (fStreamLevel==0) return 0;
 218   if (fDebugStreamer) return fDebugStreamer;
 219   TString dsName;
 220   dsName=GetName();
 221   dsName+="Debug.root";
 222   dsName.ReplaceAll(" ","");
 223   fDebugStreamer = new TTreeSRedirector(dsName.Data());
 224   return fDebugStreamer;
 225 }
 226
 227
 228
 229
 230
 231
 232 AliExternalTrackParam * AliMCTrackingTestTask::MakeTrack(const AliTrackReference* ref, TParticle*part)
 233 {
 234   //
 235   // Make track out of the track ref
 236   // part - TParticle used to determine chargr
 237   // the covariance matrix - equal 0 - starting from ideal MC position
 238   Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
 239   Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
 240   Double_t cv[21];
 241   for (Int_t i=0; i<21;i++) cv[i]=0;
 242   if (!part->GetPDG()) return 0;
 243   AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
 244   return param;
 245 }
 246
 247 Bool_t  AliMCTrackingTestTask::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
 248   //
 249   // Propagate track to point xyz using
 250   // AliTracker::PropagateTo functionality
 251   //
 252   //  param - track parameters
 253   //  xyz   - position to propagate
 254   //  mass  - particle mass
 255   //  step  - step to be used
 256   Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 257   AliTracker::PropagateTrackTo(param, radius+step, mass, step, kTRUE,0.99);
 258   AliTracker::PropagateTrackTo(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
 259   Double_t sxyz[3]={0,0,0};
 260   AliESDVertex vertex(xyz,sxyz);
 261   Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10);
 262   return isOK;
 263 }
 264
 265
 266 void  AliMCTrackingTestTask::ProcessMCInfo(){
 267   //
 268   //
 269   //
 270    //
 271   TParticle * particle= new TParticle;
 272   TClonesArray * trefs = new TClonesArray("AliTrackReference");
 273   const Double_t kPcut=0.1;
 274   //
 275   //
 276   // Process tracks
 277   //
 278   Int_t npart = fMCinfo->GetNumberOfTracks();
 279   if (npart==0) return;
 280   Double_t vertex[4]={0,0,0,0};
 281   fMCinfo->GetParticleAndTR(0, particle, trefs);
 282   if (particle){
 283     vertex[0]=particle->Vx();
 284     vertex[1]=particle->Vy();
 285     vertex[2]=particle->Vz();
 286     vertex[3]=particle->R();
 287   }
 288   //
 289   //
 290
 291   for (Int_t ipart=0;ipart<npart;ipart++){
 292     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
 293     if (status<0) continue;
 294     if (!particle) continue;
 295     if (!trefs) continue;
 296     Int_t nref = trefs->GetEntries();
 297     if (nref<5) continue;
 298     FitTrackRefs(particle,trefs);
 299     AliTrackReference * tpcIn=0;
 300     AliTrackReference * tpcOut=0;
 301     AliTrackReference * trdIn=0;
 302     AliTrackReference * trdOut=0;
 303     AliTrackReference * itsIn=0;
 304     AliTrackReference * itsOut=0;
 305     Double_t rmax=0;
 306     Double_t rmin=1000;
 307     for (Int_t iref=0;iref<nref;iref++){
 308       AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 309       if (!ref) continue;
 310
 311       Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 312
 313       if (dir<0) break; // oposite direction - looping track - return back
 314       if (ref->P()<kPcut) continue;
 315       if (ref->R()<rmax) break;
 316       //if (ref->R()<rmin)  break;
 317       //TPC
 318       if (ref->DetectorId()==AliTrackReference::kTPC){
 319         if (!tpcIn) {
 320           tpcIn  = ref;
 321         }else{
 322           if (ref->R()>tpcIn->R()) tpcOut = ref;
 323         }
 324       }
 325       //ITS
 326       if (ref->DetectorId()==AliTrackReference::kITS){
 327         if (!itsIn) {
 328           itsIn  = ref;
 329         }else{
 330           if (ref->R()>itsIn->R()) itsOut = ref;
 331         }
 332       }
 333       //TRD
 334       if (ref->DetectorId()==AliTrackReference::kTRD){
 335         if (!trdIn) {
 336           trdIn  = ref;
 337         }else{
 338           if (ref->R()>trdIn->R()) trdOut = ref;
 339         }
 340       }
 341       if (ref->R()<rmin) rmin=ref->R();
 342       if (ref->R()>rmax) rmax=ref->R();
 343     }
 344     if (tpcIn && tpcOut) {
 345       ProcessRefTracker(tpcIn,tpcOut,particle,1);
 346       ProcessRefTracker(tpcIn,tpcOut,particle,3);
 347     }
 348     if (itsIn && itsOut) ProcessRefTracker(itsIn,itsOut,particle,0);
 349     if (trdIn && trdOut) ProcessRefTracker(trdIn,trdOut,particle,2);
 350   }
 351 }
 352
 353
 354
 355 void AliMCTrackingTestTask::ProcessRefTracker(AliTrackReference* refIn,  AliTrackReference* refOut, TParticle*part,Int_t type){
 356   //
 357   // Test propagation from In to out
 358   //
 359   AliExternalTrackParam *param = 0;
 360   AliExternalTrackParam *paramMC = 0;
 361   Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
 362   Double_t mass = part->GetMass();
 363   Double_t step=1;
 364   //
 365   param=MakeTrack(refOut,part);
 366   paramMC=MakeTrack(refOut,part);
 367   if (!param) return;
 368   if (type<3) PropagateToPoint(param,xyzIn, mass, step);
 369   if (type==3) {
 370     AliTPCseed seed;
 371     seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
 372     Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
 373     seed.Rotate(alpha-seed.GetAlpha());
 374     seed.SetMass(mass);
 375     for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
 376       seed.PropagateTo(xlayer);
 377     }
 378     seed.PropagateTo(refIn->R());
 379     param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
 380   }
 381   TTreeSRedirector *pcstream = GetDebugStreamer();
 382   TVectorD gpos(3);
 383   TVectorD gmom(3);
 384   param->GetXYZ(gpos.GetMatrixArray());
 385   param->GetPxPyPz(gmom.GetMatrixArray());
 386   if (pcstream){
 387     (*pcstream)<<"MC"<<
 388       "type="<<type<<
 389       "step="<<step<<
 390       "refIn.="<<refIn<<
 391       "refOut.="<<refOut<<
 392       "p.="<<part<<
 393       "par.="<<param<<
 394       "parMC.="<<paramMC<<
 395       "gpos.="<<&gpos<<
 396       "gmom.="<<&gmom<<
 397       "\n";
 398   }
 399 }
 400
 401
 402 void  AliMCTrackingTestTask::FitTrackRefs(TParticle * part, TClonesArray * trefs){
 403   //
 404   //
 405   //
 406   //
 407   const Int_t kMinRefs=6;
 408   Int_t nrefs = trefs->GetEntries();
 409   if (nrefs<kMinRefs) return; // we should have enough references
 410   Int_t iref0 =-1;
 411   Int_t iref1 =-1;
 412
 413   for (Int_t iref=0; iref<nrefs; iref++){
 414     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 415     if (!ref) continue;
 416     Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
 417     if (dir<0) break;
 418     if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
 419     if (iref0<0) iref0 = iref;
 420     iref1 = iref;
 421   }
 422   if (iref1-iref0<kMinRefs) return;
 423   Double_t covar[14];
 424   for (Int_t icov=0; icov<14; icov++) covar[icov]=0;
 425   covar[0]=1;
 426   covar[2]=1;
 427   covar[5]=1;
 428   covar[9]=1;
 429   covar[14]=1;
 430
 431   AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
 432   AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
 433   AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
 434   AliExternalTrackParam *paramUpdate   = MakeTrack(refIn,part);
 435   paramUpdate->AddCovariance(covar);
 436   Double_t mass = part->GetMass();
 437   Double_t charge = part->GetPDG()->Charge()/3.;
 438 /*
 439   Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
 440   Float_t radiusIn= refIn->R();
 441   Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
 442   Float_t radiusOut= refOut->R();
 443 */
 444   Bool_t isOKP=kTRUE;
 445   Bool_t isOKU=kTRUE;
 446   AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
 447   for (Int_t iref = iref0; iref<=iref1; iref++){
 448     AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
 449     Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
 450     Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
 451     Double_t mag[3];
 452     field->Field(pos,mag);
 453     isOKP&=paramPropagate->Rotate(alphaC);
 454     isOKU&=paramUpdate->Rotate(alphaC);
 455     for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
 456       isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
 457       isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
 458     }
 459     isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
 460     isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
 461     Double_t clpos[2] = {0, ref->Z()};
 462     Double_t clcov[3] = { 0.005,0,0.005};
 463     isOKU&= paramUpdate->Update(clpos, clcov);
 464   }
 465   TTreeSRedirector *pcstream = GetDebugStreamer();
 466   if (pcstream){
 467     TVectorD gposU(3);
 468     TVectorD gmomU(3);
 469     TVectorD gposP(3);
 470     TVectorD gmomP(3);
 471     paramUpdate->GetXYZ(gposU.GetMatrixArray());
 472     paramUpdate->GetPxPyPz(gmomU.GetMatrixArray());
 473     paramPropagate->GetXYZ(gposP.GetMatrixArray());
 474     paramPropagate->GetPxPyPz(gmomP.GetMatrixArray());
 475
 476      (*pcstream)<<"MCupdate"<<
 477        "isOKU="<<isOKU<<
 478        "isOKP="<<isOKP<<
 479        "m="<<mass<<
 480        "q="<<charge<<
 481        "part.="<<part<<
 482        "refIn.="<<refIn<<
 483        "refOut.="<<refOut<<
 484        "pP.="<<paramPropagate<<
 485        "pU.="<<paramUpdate<<
 486        "gposU.="<<&gposU<<
 487        "gmomU.="<<&gmomU<<
 488        "gposP.="<<&gposP<<
 489        "gmomP.="<<&gmomP<<
 490        "\n";
 491    }
 492 }
 493
 494
 495
 496
 497 void AliMCTrackingTestTask::FinishTaskOutput()
 498 {
 499   //
 500   // According description in AliAnalisysTask this method is call
 501   // on the slaves before sending data
 502   //
 503   Terminate("slave");
 504   gSystem->Exec("pwd");
 505   RegisterDebugOutput();
 506
 507 }
 508
 509
 510 void AliMCTrackingTestTask::RegisterDebugOutput(){
 511   //
 512   //
 513   //
 514   //
 515   // store  - copy debug output to the destination position
 516   // currently ONLY for local copy
 517   TString dsName;
 518   dsName=GetName();
 519   dsName+="Debug.root";
 520   dsName.ReplaceAll(" ","");
 521   TString dsName2=fDebugOutputPath.Data();
 522   gSystem->MakeDirectory(dsName2.Data());
 523   dsName2+=gSystem->HostName();
 524   gSystem->MakeDirectory(dsName2.Data());
 525   dsName2+="/";
 526   dsName2+=gSystem->BaseName(gSystem->pwd());
 527   dsName2+="/";
 528   gSystem->MakeDirectory(dsName2.Data());
 529   dsName2+=dsName;
 530   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
 531   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
 532   TFile::Cp(dsName.Data(),dsName2.Data());
 533 }
 534
 535
 536 /*
 537   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 538   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
 539   AliXRDPROOFtoolkit tool;
 540   TChain * chain = tool.MakeChain("mctracking.txt","MC",0,100);
 541   chain->Lookup();
 542   //
 543   //
 544   chain->SetAlias("pdg","(p.fPdgCode)");
 545   chain->SetAlias("dPRec","(refOut.P()-par.P())/refIn.P()");
 546   chain->SetAlias("dPMC","(refOut.P()-refIn.P())/refIn->P()");
 547   chain->SetAlias("dPtRec","(refOut.Pt()-par.Pt())/refIn.Pt()");
 548   chain->SetAlias("dPtMC","(refOut.Pt()-refIn.Pt())/refIn->Pt()");
 549
 550
 551   // ITS
 552   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==0&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 553   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
 554   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
 555   gPad->SaveAs("picLoss/dPcorr_ITS_step1.gif");
 556   gPad->SaveAs("picLoss/dPcorr_ITS_step1.eps");
 557   // TPC
 558   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==1&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 559   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
 560   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
 561   gPad->SaveAs("picLoss/dPcorr_TPC_step1.gif");
 562   gPad->SaveAs("picLoss/dPcorr_TPC_step1.eps");
 563   //
 564    // TPC
 565   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 566   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
 567   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
 568   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.gif");
 569   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.eps");
 570
 571
 572   // TRD
 573   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==2&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 574   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
 575   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
 576   gPad->SaveAs("picLoss/dPcorr_TRD_step1.gif");
 577   gPad->SaveAs("picLoss/dPcorr_TRD_step1.eps");
 578
 579   //
 580   //
 581   //
 582   chain->Draw("(par.Pt()-refIn.Pt())/refIn.Pt()>>his(100,-0.02,0.02)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 583   his->SetXTitle("(P_{trec}-P_{tmc})/P_{tmc}");
 584   gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.eps");
 585   gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.gif");
 586
 587   chain->Draw("(par.P()-refIn.P())/refIn.P()>>his(100,-0.02,0.02)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
 588   his->SetXTitle("(P_{rec}-P_{mc})/P_{mc}");
 589   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.eps");
 590   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.gif");
 591 */