TPC/AliTPCcalibMaterial.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   // Load libraries
  19
  20   gSystem->Load("libANALYSIS");
  21   gSystem->Load("libTPCcalib");
  22
  23
  24   .x ~/NimStyle.C
  25   gSystem->Load("libANALYSIS");
  26   gSystem->Load("libTPCcalib");
  27
  28   // analyze results
  29
  30   TFile f("CalibObjectsTrain2.root");
  31   AliTPCcalibMaterial *calibMaterial = (AliTPCcalibMaterial *)f->Get("alignMaterial");
  32
  33
  34 */
  35
  36 #include "Riostream.h"
  37 #include "TChain.h"
  38 #include "TTree.h"
  39 #include "TH1F.h"
  40 #include "TH2F.h"
  41 #include "TH3F.h"
  42 #include "THnSparse.h"
  43 #include "TList.h"
  44 #include "TMath.h"
  45 #include "TCanvas.h"
  46 #include "TFile.h"
  47 #include "TF1.h"
  48 #include "TVectorD.h"
  49 #include "TProfile.h"
  50 #include "TGraphErrors.h"
  51 #include "TCanvas.h"
  52
  53 #include "AliTPCclusterMI.h"
  54 #include "AliTPCseed.h"
  55 #include "AliESDVertex.h"
  56 #include "AliESDEvent.h"
  57 #include "AliESDfriend.h"
  58 #include "AliESDInputHandler.h"
  59 #include "AliAnalysisManager.h"
  60
  61 #include "AliTracker.h"
  62 #include "AliMagF.h"
  63 #include "AliTPCCalROC.h"
  64
  65 #include "AliLog.h"
  66
  67 #include "AliTPCcalibMaterial.h"
  68
  69 #include "TTreeStream.h"
  70 #include "AliTPCTracklet.h"
  71 #include "TTimeStamp.h"
  72 #include "AliTPCcalibDB.h"
  73 #include "AliTPCcalibLaser.h"
  74 #include "AliDCSSensorArray.h"
  75 #include "AliDCSSensor.h"
  76
  77 #include "AliKFParticle.h"
  78 #include "AliKFVertex.h"
  79
  80 #include "AliESDTZERO.h"
  81
  82 ClassImp(AliTPCcalibMaterial)
  83
  84 AliTPCcalibMaterial::AliTPCcalibMaterial():
  85   AliTPCcalibBase("calibMaterial","calibMaterial"),
  86   fHisMaterial(0),
  87   fHisMaterialRPhi(0)
  88 {
  89
  90 }
  91
  92 AliTPCcalibMaterial::AliTPCcalibMaterial(const char * name, const char * title):
  93   AliTPCcalibBase(name,title),
  94   fHisMaterial(0),
  95   fHisMaterialRPhi(0)
  96 {
  97   //
  98   //
  99   //
 100 }
 101
 102 AliTPCcalibMaterial::~AliTPCcalibMaterial(){
 103   //
 104   // delete histograms
 105   // class is owner of all histograms
 106   //
 107   if (!fHisMaterial) return;
 108   delete fHisMaterial;
 109   delete fHisMaterialRPhi;
 110   fHisMaterial=0;
 111 }
 112
 113
 114 Long64_t AliTPCcalibMaterial::Merge(TCollection *li) {
 115   //
 116   // Merge histograms
 117   //
 118   TIterator* iter = li->MakeIterator();
 119   AliTPCcalibMaterial* cal = 0;
 120
 121   while ((cal = (AliTPCcalibMaterial*)iter->Next())) {
 122     if (!cal->InheritsFrom(AliTPCcalibMaterial::Class())) {
 123       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 124       return -1;
 125     }
 126     AliTPCcalibMaterial* calib= (AliTPCcalibMaterial*)(cal);
 127     //
 128     if (!fHisMaterial) fHisMaterial=MakeHisto();
 129     fHisMaterial->Add(calib->fHisMaterial);
 130     fHisMaterialRPhi->Add(calib->fHisMaterialRPhi);
 131   }
 132   return 0;
 133 }
 134
 135
 136
 137 void AliTPCcalibMaterial::Process(AliESDEvent *event){
 138   //
 139   //
 140   //
 141   const Int_t kMinCl=40;
 142   const Float_t kMinRatio=0.7;
 143   const Float_t kMaxS=0.05;
 144   const Float_t kMinDist=5;
 145   const Double_t kStep=1.;
 146   if (!event) return;
 147   ProcessPairs(event);
 148   return;
 149   //  TTreeSRedirector * cstream =  GetDebugStreamer();
 150   //
 151   if (!fHisMaterial){
 152     MakeHisto();
 153   }
 154
 155   //
 156   // fill histogram of track prolongations
 157   Float_t dca[2];
 158   Int_t ntracks = event->GetNumberOfTracks();
 159   for (Int_t itrack=0; itrack<ntracks; itrack++){
 160     AliESDtrack *track=event->GetTrack(itrack);
 161     if (!track) continue;
 162     if (track->GetTPCNcls()<=kMinCl) continue;
 163     if ((1.+track->GetTPCNcls())/(1.+track->GetTPCNclsF())<=kMinRatio) continue;
 164     if ((1.+track->GetTPCnclsS())/(1.+track->GetTPCNcls())>kMaxS) continue;
 165     if (!track->GetInnerParam()) continue;
 166     if (track->GetKinkIndex(0)!=0) continue;
 167     //
 168     track->GetImpactParameters(dca[0],dca[1]);
 169     if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
 170     AliExternalTrackParam param(*(track->GetInnerParam()));
 171     if (!AliTracker::PropagateTrackTo(&param,90,0.0005,10,kTRUE)) continue;
 172     Double_t x[5]={0,0,0,TMath::Sqrt(TMath::Abs(param.GetP()))*param.GetSign(),TMath::Sqrt(TMath::Abs(track->GetTPCsignal()))};
 173     //
 174     //
 175     for (Float_t radius=90; radius>0; radius-=kStep){
 176       if (!AliTracker::PropagateTrackTo(&param,radius,0.0005,kStep*0.5,kTRUE)) break;
 177       if (TMath::Abs(param.GetSnp())>0.8) break;
 178       param.GetXYZ(x);
 179       Double_t weight=1./TMath::Sqrt(1.+param.GetSnp()*param.GetSnp()+param.GetTgl()*param.GetTgl());
 180       fHisMaterial->Fill(x,weight);
 181       Double_t r = TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 182       Double_t phi = TMath::ATan2(x[1],x[0]);
 183       x[0]=r;
 184       x[1]=phi;
 185       fHisMaterialRPhi->Fill(x,weight);
 186     }
 187   }
 188 }
 189
 190 THnSparse *AliTPCcalibMaterial::MakeHisto(){
 191   //
 192   // Make track prolongation histogram
 193   //
 194   //
 195   //                    gX       gY     gz   p       dEdx
 196   Int_t    bins[5]   = {100,    100,   300,  40,   100};
 197   Double_t xmin[5]   = {-100,  -100,  -300,  -2,   5};
 198   Double_t xmax[5]   = {100,    100,   300,   2,   33};
 199   TString  axisName[5]={
 200     "gx",
 201     "gy",
 202     "gz",
 203     "p",
 204     "dedx"
 205   };
 206   TString  axisTitle[5]={
 207     "x    (cm)",
 208     "y    (cm)",
 209     "z    (cm)",
 210     "p    (GeV)",
 211     "dedx (a.u)"
 212   };
 213
 214   Int_t    binsR[5]   = {30,    360,     300,  40,   100};
 215   Double_t xminR[5]   = { 0,    -3.14,  -300,  -2,   5};
 216   Double_t xmaxR[5]   = {30,    3.14,    300,   2,   33};
 217   TString  axisNameR[5]={
 218     "r",
 219     "rphi",
 220     "z",
 221     "p",
 222     "dedx"
 223   };
 224   TString  axisTitleR[5]={
 225     "r    (cm)",
 226     "rphi    (cm)",
 227     "z    (cm)",
 228     "p    (GeV)",
 229     "dedx (a.u)"
 230   };
 231
 232   THnSparse *sparse = new THnSparseF("his_Material", "His Material", 5, bins, xmin, xmax);
 233   THnSparse *sparseR = new THnSparseF("his_MaterialRPhi", "His Material Rphi", 5, binsR, xminR, xmaxR);
 234   for (Int_t iaxis=0; iaxis<5; iaxis++){
 235     sparse->GetAxis(iaxis)->SetName(axisName[iaxis]);
 236     sparse->GetAxis(iaxis)->SetTitle(axisTitle[iaxis]);
 237     sparseR->GetAxis(iaxis)->SetName(axisNameR[iaxis]);
 238     sparseR->GetAxis(iaxis)->SetTitle(axisTitleR[iaxis]);
 239   }
 240   fHisMaterial=sparse;
 241   fHisMaterialRPhi=sparseR;
 242   return sparse;
 243 }
 244
 245
 246
 247
 248 void  AliTPCcalibMaterial::ProcessPairs(AliESDEvent *event){
 249   //
 250   // Process pairs of tracks to get a material budget map
 251   //
 252   //
 253   if (!event) return;
 254   if (event) AliKFParticle::SetField(event->GetMagneticField());  // set mean magnetic field for KF particles
 255   //
 256   // 1. Calculate total dEdx for all TPC tracks
 257   //
 258   const Int_t kMinCl=70;
 259   const Double_t kEpsilon=0.000001;
 260   const Float_t kMinRatio=0.7;
 261   const Float_t kMinDist=1.5;
 262   const Float_t kMinDistChi2=8;        //
 263   const Float_t kMaxDistZ=280;         // max distanceZ
 264   const Float_t kMaxDistR=250;          // max distanceR
 265   const Double_t kMaxChi2   =36;     // maximal chi2 to define the vertex
 266   const Double_t kMaxDistVertexSec=3;     // maximal distance to secondary vertex
 267
 268   if (!event) return;
 269   AliESDVertex *vertexSPD =  (AliESDVertex *)event->GetPrimaryVertexSPD();
 270   AliESDVertex * spdVertex   = (AliESDVertex *)event->GetPrimaryVertexSPD();
 271   AliESDVertex * trackVertex = (AliESDVertex *)event->GetPrimaryVertexTracks();
 272   AliESDVertex * tpcVertex   = (AliESDVertex *)event->GetPrimaryVertexTPC();
 273   AliESDTZERO  * tzero       = (AliESDTZERO  *)event->GetESDTZERO() ;
 274   //
 275   Double_t tpcSignalTotPrim=0;
 276   Double_t tpcSignalTotSec=0;
 277   Int_t ntracksTPC=0;
 278   Int_t nTPCPrim=0;
 279   Int_t nTPCSec=0;
 280   Int_t ntracks=event->GetNumberOfTracks();
 281   if ( ntracks<=2 ) return;
 282
 283   //
 284   Float_t dca[2]={0};
 285   Float_t cov[3]={0};
 286   Float_t dca0[2]={0};
 287   Float_t dca1[2]={0};
 288   //
 289   //1. Calculate total dEdx for primary and secondary tracks
 290   //   and count primaries and secondaries
 291   Int_t *rejectTrack = new Int_t[ntracks];
 292
 293   for (Int_t itrack=0; itrack<ntracks; itrack++){
 294     AliESDtrack *track=event->GetTrack(itrack);
 295     rejectTrack[itrack]=0;
 296     if (!track) continue;
 297     if (track->GetTPCNcls()<=kMinCl) continue;  // skip short tracks
 298     ntracksTPC++;
 299     if ((1.+track->GetTPCNcls())/(1.+track->GetTPCNclsF())<=kMinRatio) continue;
 300     if (!track->GetInnerParam()) continue;  // skip not TPC tracks
 301     if (track->GetKinkIndex(0)!=0)  {rejectTrack[itrack]+=16;continue;} // skip kinks
 302     track->GetImpactParameters(dca[0],dca[1]);
 303     if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
 304     // remove too dip secondaries
 305     if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist){
 306       tpcSignalTotPrim+=track->GetTPCsignal();
 307       nTPCPrim++;
 308     }else{
 309       tpcSignalTotSec+=track->GetTPCsignal();
 310       nTPCSec++;
 311     };
 312     if (cov[0]>kEpsilon &&TMath::Abs(dca[0])>kEpsilon &&TMath::Sqrt(dca[0]*dca[0]/(TMath::Abs(cov[0])))<kMinDistChi2) rejectTrack[itrack]+=1;  // primary
 313     if (cov[0]>kEpsilon &&TMath::Abs(dca[0])>kEpsilon &&TMath::Abs(dca[0])<kMinDist)  rejectTrack[itrack]+=1;  // primary
 314     if (track->GetTPCsignal()<40) rejectTrack[itrack]+=16;
 315     //
 316     if (CheckLooper(itrack, event))   rejectTrack[itrack]+=2;   // looper
 317     if (CheckV0(itrack,event))       rejectTrack[itrack]+=4;
 318   }
 319
 320
 321   //
 322   // 2. Find secondary vertices - double loop
 323   //
 324   for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
 325     AliESDtrack *track0=event->GetTrack(itrack0);
 326     if (!track0) continue;
 327     if (track0->GetTPCNcls()<=kMinCl) continue;  // skip short tracks
 328     if ((1.+track0->GetTPCNcls())/(1.+track0->GetTPCNclsF())<=kMinRatio) continue;
 329     if (!track0->GetInnerParam()) continue;  // skip not TPC tracks
 330     if (track0->GetKinkIndex(0)>0) continue; // skip kinks
 331     if (rejectTrack[itrack0]) continue;   // skip
 332     track0->GetImpactParameters(dca[0],dca[1]);
 333     track0->GetImpactParameters(dca0[0],dca0[1]);
 334     if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
 335     // remove too dip secondaries
 336
 337     AliKFParticle part0(*track0,211);  //assuming pion mass
 338     if (track0->Charge()*part0.Q()<0) part0.Q()*=-1;  // change sign if opposite
 339     //
 340     for (Int_t itrack1=itrack0+1; itrack1<ntracks; itrack1++){
 341       AliESDtrack *track1=event->GetTrack(itrack1);
 342       if (!track1) continue;
 343       if (rejectTrack[itrack1]) continue;   // skip
 344       if (track1->GetTPCNcls()<=kMinCl) continue;  // skip short tracks
 345       if ((1.+track1->GetTPCNcls())/(1.+track1->GetTPCNclsF())<=kMinRatio) continue;
 346       if (!track1->GetInnerParam()) continue;  // skip not TPC tracks
 347       if (track1->GetKinkIndex(0)!=0) continue; // skip kinks
 348       track1->GetImpactParameters(dca1[0],dca1[1]);
 349       track1->GetImpactParameters(dca[0],dca[1]);
 350       if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
 351       if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
 352       AliKFParticle part1(*track1,211); // assuming pion mass
 353       if (track1->Charge()*part1.Q()<0) part1.Q()*=-1;  // change sign if opposite
 354
 355       //
 356       //
 357       AliKFVertex vertex;
 358       vertex+=part0;
 359       vertex+=part1;
 360       if ((vertex.GetChi2()/vertex.GetNDF())> kMaxChi2) continue;
 361       if (TMath::Abs(vertex.GetX())>kMaxDistR) continue;
 362       if (TMath::Abs(vertex.GetY())>kMaxDistR) continue;
 363       if (TMath::Abs(vertex.GetZ())>kMaxDistZ) continue;
 364       Double_t errX2=vertex.GetErrX();
 365       Double_t errY2=vertex.GetErrY();
 366       Double_t errZ2=vertex.GetErrZ();
 367       //
 368       Double_t err3D=TMath::Sqrt(errX2*errX2+errY2*errY2+errZ2*errZ2/10.);
 369       Double_t err2D=TMath::Sqrt(errX2*errX2+errY2*errY2);
 370       if (err3D>kMaxDistVertexSec) continue;
 371       if (err3D*TMath::Sqrt(vertex.GetChi2()+0.00001)>kMaxDistVertexSec) continue;
 372
 373       Double_t dvertex=0;
 374       dvertex += (vertexSPD->GetX()-vertex.GetX())*(vertexSPD->GetX()-vertex.GetX());
 375       dvertex += (vertexSPD->GetY()-vertex.GetY())*(vertexSPD->GetY()-vertex.GetY());
 376       dvertex += (vertexSPD->GetZ()-vertex.GetZ())*(vertexSPD->GetZ()-vertex.GetZ());
 377       dvertex=TMath::Sqrt(dvertex+0.00000001);
 378       if (err3D>0.2*dvertex) continue;
 379       if (err3D*TMath::Sqrt(vertex.GetChi2()+0.000001)>0.1*dvertex) continue;
 380       Double_t radius = TMath::Sqrt((vertex.GetX()*vertex.GetX()+vertex.GetY()*vertex.GetY()));
 381       //
 382       AliKFVertex vertex2;
 383       vertex2+=part0;
 384       vertex2+=part1;
 385       //
 386
 387       for (Int_t itrack2=0; itrack2<ntracks; itrack2++){
 388         if (itrack2==itrack0) continue;
 389         if (itrack2==itrack1) continue;
 390         if (rejectTrack[itrack2]) continue;   // skip
 391         AliESDtrack *track2=event->GetTrack(itrack2);
 392         if (!track2) continue;
 393         if (track2->GetTPCNcls()<=kMinCl) continue;  // skip short tracks
 394         if ((1.+track2->GetTPCNcls())/(1.+track2->GetTPCNclsF())<=kMinRatio) continue;
 395         if (!track2->GetInnerParam()) continue;  // skip not TPC tracks
 396         if (track2->GetKinkIndex(0)>0) continue; // skip kinks
 397         track2->GetImpactParameters(dca[0],dca[1]);
 398         if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
 399         if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
 400         if (TMath::Abs(track2->GetD(vertex.GetX(), vertex.GetY(),event->GetMagneticField()))>kMaxDistVertexSec) continue;
 401         Double_t vtxx[3]={vertex2.GetX(),vertex2.GetY(),vertex2.GetZ()};
 402         Double_t svtxx[3]={vertex.GetErrX(),vertex.GetErrY(),vertex.GetErrZ()};
 403         AliESDVertex vtx(vtxx,svtxx);
 404         AliExternalTrackParam param=*track2;
 405         Double_t delta[2]={0,0};
 406         if (!param.PropagateToDCA(&vtx,event->GetMagneticField(),kMaxDistVertexSec,delta)) continue;
 407         if (TMath::Abs(delta[0])>kMaxDistVertexSec) continue;
 408         if (TMath::Abs(delta[1])>kMaxDistVertexSec) continue;
 409         if (TMath::Abs(delta[0])>6.*TMath::Sqrt(param.GetSigmaY2()+vertex2.GetErrY()*vertex2.GetErrY())+0.1) continue;
 410         if (TMath::Abs(delta[1])>6.*TMath::Sqrt(param.GetSigmaZ2()+vertex2.GetErrZ()*vertex2.GetErrZ())+0.5) continue;
 411         //
 412         AliKFParticle part2(param,211); // assuming pion mass
 413         if (track2->Charge()*part2.Q()<0) part2.Q()*=-1;  // change sign if opposite
 414         vertex2+=part2;
 415         rejectTrack[itrack0]+=10;  // do noit reuse the track
 416         rejectTrack[itrack1]+=10;  // do not reuse the track
 417         rejectTrack[itrack2]+=10;
 418       }
 419
 420
 421       TTreeSRedirector *pcstream = GetDebugStreamer();
 422       if (pcstream){
 423         //
 424         //
 425         //
 426         Float_t dedx0= track0->GetTPCsignal();
 427         Float_t dedx1= track1->GetTPCsignal();
 428         AliExternalTrackParam * p0= (AliExternalTrackParam *)track0->GetInnerParam();
 429         AliExternalTrackParam * p1= (AliExternalTrackParam *)track1->GetInnerParam();
 430         Double_t errX=vertex2.GetErrX();
 431         Double_t errY=vertex2.GetErrY();
 432         Double_t errZ=vertex2.GetErrZ();
 433         Double_t vx = vertex2.GetX();
 434         Double_t vy = vertex2.GetY();
 435         Double_t vz = vertex2.GetZ();
 436         (*pcstream)<<"mapTPC"<<
 437           "run="<<fRun<<                     // run
 438           "event="<<fEvent<<          //  event number
 439           "time="<<fTime<<         // timeStamp
 440           "trigger="<<fTrigger<<      //  trigger
 441           "mag="<<fMagF<<             //  magnetic field
 442           //
 443           "spdV.="<<spdVertex<<            // spd vertex
 444           "trackV.="<<trackVertex<<        // track vertex
 445           "tpcV.="<<tpcVertex<<            // track vertex
 446           "tzero.="<<tzero<<               // tzero info
 447           //
 448           "ntracks="<<ntracks<<
 449           "ntracksTPC="<<ntracksTPC<<
 450           "nPrim="<<nTPCPrim<<              // number of primaries
 451           "nSec="<<nTPCSec<<                // number of secondaries
 452           "sigPrim="<<tpcSignalTotPrim<<    // total dEdx in primaries
 453           "sigSec="<<tpcSignalTotSec<<      // total dEdx in secondaries
 454           "dedx0="<<dedx0<<                 // dedx part 0
 455           "dedx1="<<dedx1<<                 // dedx part 1
 456           "p0.="<<p0<<                      // part 0
 457           "p1.="<<p1<<                       //part 1
 458           "v.="<<&vertex<<                  // KF vertex
 459           "v2.="<<&vertex2<<                // KF vertex all tracks
 460           "z0="<<dca0[1]<<
 461           "z1="<<dca1[1]<<
 462           "rphi0="<<dca0[0]<<
 463           "rphi1="<<dca1[0]<<
 464           "dvertex="<<dvertex<<
 465           "radius="<<radius<<
 466           "vx="<<vx<<
 467           "vy="<<vy<<
 468           "vz="<<vz<<
 469           "errX="<<errX<<
 470           "errY="<<errY<<
 471           "errZ="<<errZ<<
 472           "err2D="<<err2D<<
 473           "err3D="<<err3D<<
 474           "\n";
 475         //
 476         if (vertex2.GetNDF()>2){
 477           (*pcstream)<<"mapVertex"<<
 478             "run="<<fRun<<                     // run
 479             "event="<<fEvent<<          //  event number
 480             "time="<<fTime<<         // timeStamp
 481             "trigger="<<fTrigger<<      //  trigger
 482             "mag="<<fMagF<<             //  magnetic field
 483             //
 484             "spdV.="<<spdVertex<<            // spd vertex
 485             "trackV.="<<trackVertex<<        // track vertex
 486             "tpcV.="<<tpcVertex<<            // track vertex
 487             "tzero.="<<tzero<<               // tzero info
 488             //
 489             //
 490             "ntracks="<<ntracks<<
 491             "ntracksTPC="<<ntracksTPC<<
 492             "nPrim="<<nTPCPrim<<              // number of primaries
 493             "nSec="<<nTPCSec<<                // number of secondaries
 494             "sigPrim="<<tpcSignalTotPrim<<    // total dEdx in primaries
 495             "sigSec="<<tpcSignalTotSec<<      // total dEdx in secondaries
 496             "dedx0="<<dedx0<<                 // dedx part 0
 497             "dedx1="<<dedx1<<                 // dedx part 1
 498             "p0.="<<p0<<                      // part 0
 499             "p1.="<<p1<<                       //part 1
 500             "v.="<<&vertex<<                  // KF vertex
 501             "v2.="<<&vertex2<<                // KF vertex all tracks
 502             "z0="<<dca0[1]<<
 503             "z1="<<dca1[1]<<
 504             "rphi0="<<dca0[0]<<
 505             "rphi1="<<dca1[0]<<
 506             "radius="<<radius<<
 507             "vx="<<vx<<
 508             "vy="<<vy<<
 509             "vz="<<vz<<
 510             "errX="<<errX<<
 511             "errY="<<errY<<
 512             "errZ="<<errZ<<
 513             "err2D="<<err2D<<
 514             "err3D="<<err3D<<
 515             "dvertex="<<dvertex<<
 516             "\n";
 517         }
 518       }
 519     }
 520   }
 521   TTreeSRedirector *pcstream = GetDebugStreamer();
 522   if (pcstream){
 523     (*pcstream)<<"statTPC"<<
 524       "run="<<fRun<<                     // run
 525       "time="<<fTime<<         // timeStamp
 526       "trigger="<<fTrigger<<      //  trigger
 527       "mag="<<fMagF<<             //  magnetic field
 528       "ntracks="<<ntracks<<
 529       "ntracksTPC="<<ntracksTPC<<
 530       //
 531       "nPrim="<<nTPCPrim<<              // number of primaries
 532       "nSec="<<nTPCSec<<                // number of secondaries
 533       "sigPrim="<<tpcSignalTotPrim<<    // total dEdx in primaries
 534       "sigSec="<<tpcSignalTotSec<<      // total dEdx in secondaries
 535       //
 536       "spdV.="<<spdVertex<<            // spd vertex
 537       "trackV.="<<trackVertex<<        // track vertex
 538       "tpcV.="<<tpcVertex<<            // track vertex
 539       "tzero.="<<tzero<<               // tzero info
 540       "\n";
 541   }
 542   delete [] rejectTrack;
 543 }
 544
 545
 546 Bool_t AliTPCcalibMaterial::CheckLooper(Int_t index, AliESDEvent *event){
 547   //
 548   // check if given track is looper candidate
 549   // if looper return kTRUE
 550   //
 551   Int_t ntracks=event->GetNumberOfTracks();
 552   Int_t index1=-1;
 553   const Double_t ktglCut=0.03;
 554   const Double_t kalphaCut=0.4;
 555   static Int_t counter=0;
 556   //
 557   AliESDtrack * track0 = event->GetTrack(index);
 558   AliESDtrack * track1P = 0;
 559   for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
 560     if (itrack1==index) continue;
 561     AliESDtrack *track1=event->GetTrack(itrack1);
 562     if (!track1) continue;
 563     if (TMath::Abs(TMath::Abs(track1->GetTgl())-TMath::Abs(track0->GetTgl()))>ktglCut) continue;
 564     if (TMath::Abs(TMath::Abs(track1->GetAlpha())-TMath::Abs(track0->GetAlpha()))>kalphaCut) continue;
 565     index1=index;
 566     track1P=track1;
 567   }
 568   if (index1>=0){
 569     TTreeSRedirector *pcstream = GetDebugStreamer();
 570     if (pcstream &&counter<100000){
 571       counter++;
 572       AliExternalTrackParam p0(*track0);
 573       AliExternalTrackParam p1(*track1P);
 574       (*pcstream)<<"checkLooper"<<
 575         "p0.="<<&p0<<
 576         "p1.="<<&p1<<
 577         "\n";
 578     }
 579     return kTRUE;
 580   }
 581   return kFALSE;
 582 }
 583
 584 Bool_t AliTPCcalibMaterial::CheckV0(Int_t index, AliESDEvent *event){
 585   //
 586   // check if given track is V0 candidata
 587   // if looper return kTRUE
 588   //
 589   return kFALSE;
 590   Int_t ntracks=event->GetNumberOfTracks();
 591   Int_t index1=-1;
 592   const Double_t kSigmaMass=0.001;
 593   const Int_t kChi2Cut=10;
 594   static Int_t counter=0;
 595   //
 596   AliESDtrack * track0 = event->GetTrack(index);
 597   AliExternalTrackParam pL(*track0);
 598   AliKFParticle part0El(*track0, 11);  //assuming  mass e
 599   AliKFParticle part0Pi(*track0, 211);  //assuming  mass pi0
 600   AliKFParticle part0P(*track0, 2212);  //assuming  mass proton
 601   if (track0->Charge()*part0El.Q()<0) {
 602     part0El.Q()*=-1;  // change sign if opposite
 603     part0Pi.Q()*=-1;  // change sign if opposite
 604     part0P.Q()*=-1;   // change sign if opposite
 605   }
 606   Bool_t isGamma=0;
 607   Bool_t isK0=0;
 608
 609   for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
 610     if (itrack1==index) continue;
 611     AliESDtrack *track1=event->GetTrack(itrack1);
 612     if (!track1) continue;
 613     if (track1->Charge()*track0->Charge()>0) continue;
 614     AliKFParticle part1El(*track1, 11);  //assuming  mass e
 615     AliKFParticle part1Pi(*track1, 211);  //assuming  mass e
 616     AliKFParticle part1P(*track1, 2212);  //assuming  mass e
 617     if (track1->Charge()*part1El.Q()<0) {
 618       part1El.Q()*=-1;  // change sign if opposite
 619       part1Pi.Q()*=-1;  // change sign if opposite
 620       part1P.Q()*=-1;   // change sign if opposite
 621     }
 622     //
 623     AliKFVertex vertexG;  // gamma conversion candidate
 624     vertexG+=part0El;
 625     vertexG+=part1El;
 626     AliKFVertex vertexGC;  // gamma conversion candidate
 627     vertexGC+=part0El;
 628     vertexGC+=part1El;
 629     vertexGC.SetMassConstraint(0,kSigmaMass);
 630     AliKFVertex vertexK0;  // gamma conversion candidate
 631     vertexK0+=part0Pi;
 632     vertexK0+=part1Pi;
 633     AliKFVertex vertexK0C;  // gamma conversion candidate
 634     vertexK0C+=part0Pi;
 635     vertexK0C+=part1Pi;
 636     vertexK0C.SetMassConstraint(0.497614,kSigmaMass);
 637     if (vertexGC.GetChi2()<kChi2Cut && vertexG.GetMass()<0.06)      isGamma=kTRUE;
 638     if (vertexK0C.GetChi2()<kChi2Cut&&TMath::Abs(vertexK0.GetMass()-0.5)<0.06)  isK0=kTRUE;
 639     if (isGamma||isK0) {
 640       index1=index;
 641       TTreeSRedirector *pcstream = GetDebugStreamer();
 642       if (pcstream&&counter<2000){
 643         counter++;
 644         AliExternalTrackParam p0(*track0);
 645         AliExternalTrackParam p1(*track1);
 646         (*pcstream)<<"checkV0"<<
 647           "p0.="<<&p0<<  //particle 0
 648           "p1.="<<&p1<<  //particle 1
 649           "isGamma="<<isGamma<<  // is gamma candidate
 650           "isK0s="<<isK0<<      // is k0 candidate
 651           "vG.="<<&vertexG<<     // is gamma candidate
 652           "vGC.="<<&vertexGC<<   // is gamma candidate
 653           "vK.="<<&vertexK0<<     // is K0s candidate
 654           "vKC.="<<&vertexK0C<<   // is K0s candidate
 655           "\n";
 656       }
 657       break;
 658     }
 659   }
 660   if (index1>0) return kTRUE;
 661   return kFALSE;
 662 }
 663
 664 /*
 665   //AliXRDPROOFtoolkit::FilterList("mater.list","* mapTPC",1);
 666   AliXRDPROOFtoolkit toolkit;
 667   TChain *chain   = toolkit.MakeChainRandom("mater.list.Good","mapVertex",0,4000);
 668   TChain *chainTPC   = toolkit.MakeChainRandom("mater.list.Good","mapTPC",0,50000);
 669
 670
 671   TCut cutErr="sqrt(v.fChi2)*err3D<1.0";
 672   TCut cutOccu="sqrt(v.fChi2*(errX^2+errY^2))/min(sqrt(v.fP[0]^2+v.fP[1]^2+v.fP[2]^2/20),20)<0.2&&sqrt((errX^2+errY^2))/min(sqrt(v.fP[0]^2+v.fP[1]^2+v.fP[2]^2/20),20)<0.3&&sqrt(v.fChi2)*err3D<1.5";
 673   //
 674   chainTPC->Draw(">>listTPC",cutOccu,"entryList");
 675   TEntryList *elistTPC = (TEntryList*)gDirectory->Get("listTPC");
 676   chainTPC->SetEntryList(elistTPC);
 677
 678   //
 679   //
 680
 681   chainTPC.Draw("v.fP[1]:v.fP[0]","sqrt(v.fP[0]^2+v.fP[1]^2)<100","",10000000);
 682
 683   TCut cutITS="abs(v.fP[2])-10<sqrt(v.fP[0]^2+v.fP[1]^2)";
 684   chainTPC.Draw("v.fP[1]:v.fP[0]",cutITS+"(sqrt(v.fP[0]^2+v.fP[1]^2)<50)&&err2D<1&&err3D/dvertex<0.05","",500000);
 685
 686
 687   chainTPC.Draw("v.fP[1]:v.fP[0]>>his(400,-100,100,400,-100,100)","err2D<1&&err2D/dvertex<0.2","colz",50000);
 688
 689
 690   chainTPC.Draw("radius:vz>>his(400,-100,100,400,0,100)","err2D<1&&err2D/dvertex<0.05","colz",5000);
 691
 692   TTree * tree = chain->CloneTree();
 693   TFile f("material.root","recreate");
 694   tree->Write();
 695   f.Close();
 696
 697
 698
 699   TCut cutChi2="v.fChi2<2&&sqrt(errX^2+errY^2)<2.";
 700   TFile f("material.root")
 701   TTree * tree = (TTree*)f.Get("mapVertex");
 702
 703   tree->Draw("v.fChi2",cutChi2,"",1000)
 704
 705
 706
 707
 708 */
 709
 710