TPC/Upgrade/AliToyMCDrawer.cxx

   1 #include <iostream>
   2 #include <TROOT.h>
   3 #include <TChain.h>
   4 #include <TMath.h>
   5 #include <TTree.h>
   6 #include <TFile.h>
   7 #include <TGraph2D.h>
   8 #include <TH3F.h>
   9 #include <TH2F.h>
  10 #include <TClonesArray.h>
  11 #include <TPolyLine3D.h>
  12 #include <TPolyMarker3D.h>
  13 #include <TAxis.h>
  14 #include <TLegend.h>
  15 #include <TPad.h>
  16 #include <TGeoGlobalMagField.h>
  17
  18 #include <AliMagF.h>
  19 #include <AliCDBManager.h>
  20 #include <AliTPCParam.h>
  21 #include <AliGeomManager.h>
  22 #include <AliTPCcalibDB.h>
  23 #include <AliTrackPointArray.h>
  24 #include <AliCluster.h>
  25 #include <AliLog.h>
  26 #include <AliTPCLaserTrack.h>
  27
  28
  29 #include "AliToyMCDrawer.h"
  30 #include "AliToyMCEvent.h"
  31 #include "AliToyMCTrack.h"
  32
  33 // Visualization class. To use
  34
  35 //  AliToyMCDrawer* draw = new AliToyMCDrawer()
  36 //  draw->SetFileName("path/to/toyMC.root")
  37
  38 //  draw->FillEventArray(Int_t centerEventNumber)
  39 //         or
  40 //  draw->FillEventArray(Double_t time)
  41 //    to display with a certain event in the center or at a certain time
  42
  43 //  draw->DrawEvents(Bool_t both, Bool_t before)
  44 //    where "both" will display events before and after the middle event and
  45 //    before will show also events before (after) the middle event if true (false)
  46 //    when "both" is false
  47
  48
  49 ClassImp(AliToyMCDrawer);
  50
  51 AliToyMCDrawer::AliToyMCDrawer()
  52   : TObject()
  53   ,fInputTree(0x0)
  54   ,fInFile(0x0)
  55   ,fFileName()
  56   ,fEvent(0x0)
  57   ,fEventArray(0x0)
  58   ,fDispHist(0x0)
  59   ,fCenterTime(-1.)
  60   ,fDriftVel(-1.)
  61   ,fTPCParam(0x0)
  62   ,fMaxZ0(0.)
  63   ,fIFCRadius(83.5)
  64   ,fOFCRadius(254.5)
  65   ,fTimeRange(0)
  66   ,fRoc(AliTPCROC::Instance())
  67   ,fPoints(0x0)
  68   ,fDistPoints(0x0)
  69   ,fProjectionType("XYT")
  70   ,fTimeZmin(0.)
  71   ,fTimeZmax(0.)
  72   ,fGlobalXmin(0.)
  73   ,fGlobalXmax(0.)
  74   ,fGlobalYmin(0.)
  75   ,fGlobalYmax(0.)
  76   {
  77    fEventArray = new TClonesArray("AliToyMCEvent");
  78
  79    fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
  80    fTPCParam->ReadGeoMatrices();
  81    fDriftVel = fTPCParam->GetDriftV();
  82    fMaxZ0    =fTPCParam->GetZLength();
  83    fTimeRange = 2*fMaxZ0/fDriftVel;
  84
  85    fPoints = new TClonesArray("TPolyMarker3D");
  86    fDistPoints = new TClonesArray("TPolyMarker3D");
  87  }
  88 //________________________________________________________________
  89 AliToyMCDrawer::AliToyMCDrawer(const AliToyMCDrawer &drawer)
  90   : TObject(drawer)
  91   ,fInputTree(drawer.fInputTree)
  92   ,fInFile(drawer.fInFile)
  93   ,fFileName(drawer.fFileName)
  94   ,fEvent(drawer.fEvent)
  95   ,fEventArray(drawer.fEventArray)
  96   ,fDispHist(drawer.fDispHist)
  97   , fCenterTime(drawer.fCenterTime)
  98   ,fDriftVel(drawer.fDriftVel)
  99   ,fTPCParam(drawer.fTPCParam)
 100   ,fMaxZ0(drawer.fMaxZ0)
 101   ,fIFCRadius(drawer.fIFCRadius)
 102   ,fOFCRadius(drawer.fOFCRadius)
 103   ,fTimeRange(drawer.fTimeRange)
 104   ,fRoc(drawer.fRoc)
 105   ,fPoints(drawer.fPoints)
 106   ,fDistPoints(drawer.fDistPoints)
 107   ,fProjectionType("XYT")
 108   ,fTimeZmin(0.)
 109   ,fTimeZmax(0.)
 110   ,fGlobalXmin(0.)
 111   ,fGlobalXmax(0.)
 112   ,fGlobalYmin(0.)
 113   ,fGlobalYmax(0.)
 114   {
 115   //
 116 }
 117 //_____________________________________________________
 118 AliToyMCDrawer& AliToyMCDrawer::operator = (const AliToyMCDrawer &drawer)
 119 {
 120   //assignment operator
 121   if (&drawer == this) return *this;
 122   new (this) AliToyMCDrawer(drawer);
 123
 124   return *this;
 125 }
 126
 127 //________________________________________________________________
 128 AliToyMCDrawer::~AliToyMCDrawer()
 129 {
 130   //destructor
 131   delete fEvent;
 132   delete fEventArray;
 133   delete fPoints;
 134   delete fDistPoints;
 135   delete fDispHist;
 136 }
 137 //________________________________________________________________
 138 Int_t AliToyMCDrawer::FillEventArray(Double_t snapShotTime)
 139 {
 140   if(fFileName.IsNull()) {
 141     std::cout << "no input file provided, using default (toyMC.root)" << std::endl;
 142     fFileName = "toyMC.root";
 143   }
 144
 145   fInFile = new TFile(fFileName.Data(),"read");
 146   gROOT->cd();
 147   fInputTree = dynamic_cast<TTree*> (fInFile->Get("toyMCtree"));
 148   fInputTree->SetBranchAddress("event",&fEvent);
 149   fEventArray->Clear();
 150
 151   Int_t leftSearchIndex  = Int_t(snapShotTime/2e-5);
 152   if(leftSearchIndex>fInputTree->GetEntries()) {
 153     leftSearchIndex = fInputTree->GetEntries()-1;
 154     fInputTree->GetEvent(leftSearchIndex);
 155     snapShotTime = fEvent->GetT0();
 156     std::cout << "input time too large, setting time to time of last event" << std::endl;
 157   }
 158   if(leftSearchIndex<0) leftSearchIndex = 0;
 159
 160   //fCenterTime = snapShotTime;
 161
 162   fInputTree->GetEvent(leftSearchIndex);
 163   Double_t leftTime = fEvent ->GetT0();
 164   Int_t firstEventIndex;
 165   //if fEvent exactly at snapShotTime, check if there are more events at the same time but with lower indices
 166   if( TMath::Abs(leftTime - snapShotTime)<5e-9)
 167     {
 168       //  printf("close from start, lefttime = %2.2f\n",leftTime);
 169       firstEventIndex = leftSearchIndex+1;
 170       while( TMath::Abs(leftTime- snapShotTime) < 5e-9){
 171         firstEventIndex--;
 172         //      printf("close from start, lefttime = %2.2f\n",leftTime);
 173         if(fInputTree->GetEvent(firstEventIndex-1))
 174           {
 175             leftTime=fEvent->GetT0();
 176             //     printf("close from start, lefttime = %2.2f\n",leftTime);
 177           }
 178         else break;
 179       }
 180
 181
 182     }
 183   else
 184     {
 185       Int_t rightSearchIndex = leftSearchIndex;
 186       Double_t rightTime;
 187
 188       Int_t direction = 1; //go right
 189       if(leftTime > snapShotTime) {
 190         direction = -1;
 191       }
 192       rightSearchIndex+= direction;
 193       if(!fInputTree->GetEvent(rightSearchIndex)) return FillEventArray(leftSearchIndex);
 194       rightTime = fEvent->GetT0();
 195       // std::cout << "bef search " << leftTime << " " << rightTime<< " " << " " << snapShotTime << " "<< (leftTime-snapShotTime)*(rightTime-snapShotTime)  << std::endl;
 196       // Int_t b;
 197
 198       while ( (leftTime-snapShotTime)*(rightTime-snapShotTime)>0   ){
 199         //              printf("searching, leftime %f, righttim %f\n",leftTime,rightTime);
 200         rightSearchIndex += direction;
 201         leftSearchIndex +=direction;
 202         fInputTree->GetEvent(leftSearchIndex);
 203         leftTime = fEvent->GetT0();
 204         if(!fInputTree->GetEvent(rightSearchIndex)) {
 205           rightSearchIndex-=direction;
 206           break;
 207
 208         }
 209         rightTime = fEvent->GetT0();
 210                 printf("searching, leftime %f, righttim %f\n",leftTime,rightTime);
 211       }
 212       if (direction==-1) rightSearchIndex = leftSearchIndex;
 213       firstEventIndex = rightSearchIndex;
 214
 215     }
 216   // fInputTree->GetEvent(firstEventIndex);
 217   //std::cout <<"first event after sn time: " << fEvent->GetT0() << std::endl;
 218   return FillEventArray(firstEventIndex, snapShotTime);
 219
 220 }
 221 //________________________________________________________________
 222 Int_t AliToyMCDrawer::FillEventArray(Int_t middleEventNbr, Double_t snapShotTime)
 223 {
 224   if(fFileName.IsNull()) {
 225     std::cout << "no input file provided, using default (toyMC.root)" << std::endl;
 226     fFileName = "toyMC.root";
 227   }
 228   if(!fInFile) fInFile = new TFile(fFileName.Data(),"read");
 229   gROOT->cd();
 230   if(!fInputTree)
 231     {
 232       fInputTree = dynamic_cast<TTree*> (fInFile->Get("toyMCtree"));
 233       //   fInputTree->Add(fileName);
 234       fInputTree->SetBranchAddress("event",&fEvent);
 235     }
 236   Double_t centerEventTime;
 237   if(fInputTree->GetEvent(middleEventNbr))
 238     {
 239       centerEventTime = fEvent->GetT0();
 240
 241     }
 242   else return 0;
 243
 244   if(snapShotTime<0) fCenterTime = centerEventTime;
 245   else fCenterTime = snapShotTime;
 246
 247   fEventArray->Clear();
 248   if(fInputTree->GetEvent(middleEventNbr)){
 249     new((*fEventArray)[fEventArray->GetEntriesFast()]) AliToyMCEvent(*fEvent);
 250
 251     //add events after middle event
 252     Int_t eventIndex = middleEventNbr + 1;
 253     Double_t currentTime = centerEventTime;
 254     if(fInputTree->GetEvent(eventIndex)) {
 255       currentTime = fEvent->GetT0();
 256       while(currentTime - centerEventTime < fTimeRange/2) {
 257         new((*fEventArray)[fEventArray->GetEntriesFast()]) AliToyMCEvent(*fEvent);
 258         eventIndex++;
 259         if(!fInputTree->GetEvent(eventIndex)) break;
 260         currentTime = fEvent->GetT0();
 261       }
 262     }
 263     //add events before middle event
 264     eventIndex = middleEventNbr - 1;
 265     if(fInputTree->GetEvent(eventIndex)){
 266       currentTime = fEvent->GetT0();
 267       while(centerEventTime-currentTime < fTimeRange/2) {
 268
 269         new((*fEventArray)[fEventArray->GetEntriesFast()]) AliToyMCEvent(*fEvent);
 270         eventIndex--;
 271         if(!fInputTree->GetEvent(eventIndex)) break;
 272         currentTime = fEvent->GetT0();
 273       }
 274     }
 275
 276     std::cout << " end fill" << std::endl;
 277
 278     return 1;
 279   }
 280   else {
 281     printf("Selected event number (%d) out of range!\n",middleEventNbr);
 282      return 0;
 283   }
 284
 285 }
 286 //________________________________________________________________
 287 void AliToyMCDrawer::DrawEvents(Bool_t both, Bool_t before)
 288 {
 289
 290   fPoints->Clear();
 291   fDistPoints->Clear();
 292   DrawGeometry();
 293
 294
 295    Double_t phiCut = 1.;
 296    TLegend *leg = new TLegend(0.1,0.75,0.6,1,Form("Snapshot of TPC at %f micros after sim start, freq. = , bunchcr.rate = .",1000000*fCenterTime));
 297    leg->AddEntry((TObject*)0,Form("%2.2f<#phi<#pi-%2.2f && #pi+%2.2f<#phi<2#pi-%2.2f",phiCut,phiCut,phiCut,phiCut),"");
 298
 299
 300
 301
 302   for(Int_t iEvent = 0; iEvent<fEventArray->GetEntriesFast(); iEvent++){
 303
 304
 305      AliToyMCEvent *currentEvent = static_cast<AliToyMCEvent*> (fEventArray->At(iEvent));
 306      Double_t currentEventTime = currentEvent->GetT0();
 307
 308      if((1000000*currentEventTime)>((1000000*fCenterTime)+0.05   ))  {
 309        printf("larger iEvent: %d, current %f, center %f, ntracks: %d\n",iEvent, currentEventTime,fCenterTime, currentEvent->GetNumberOfTracks());
 310        printf("after\n");
 311        if ( !(both || !before)) continue;
 312
 313      }
 314      if( currentEventTime<=fCenterTime)// && !(currentEventTime==fCenterTime &&both   ))
 315        {
 316          printf("smaller iEvent: %d, current %f, center %f, ntracks: %d\n",iEvent, currentEventTime,fCenterTime, currentEvent->GetNumberOfTracks());
 317          printf("before\n");
 318          if ( !(both || before)) continue;
 319        }
 320
 321      //std:: cout << "iev " << iEvent << " " << color << std::endl;
 322      Double_t currentEventTimeRelToCentral = fCenterTime - currentEventTime;
 323
 324      Int_t color = ((currentEvent->GetEventNumber())%10);
 325      if(color == 0||color ==1) {
 326        color==0?color=2:color=3;
 327        // if(iEvent==0) color = 8;
 328        // if(iEvent==1) color = 9;
 329
 330      }
 331
 332       char bef[] = "before snapshot time";
 333       char aft[] = "after snapshot time";
 334       char *when;
 335       Int_t sign = 1;
 336       if( 1000000*(currentEventTime-fCenterTime) > 0.5) {when = aft;
 337         sign = -1;
 338           }
 339       else when = bef;
 340       TGraph *temp = new TGraph();
 341       temp->SetName(Form("temp%d",currentEvent->GetEventNumber()));
 342        temp->SetLineColor(color);
 343        temp->SetLineWidth(10);
 344
 345        leg->AddEntry(temp,Form(" %2.2f microseconds %s (global time %2.2f) at z_vertex = %2.2f, centr = ., eventNr: %d", 1000000*sign*currentEventTimeRelToCentral, when, 1000000*currentEvent->GetT0(),  currentEvent->GetZ() ,currentEvent->GetEventNumber() ),"l");
 346
 347
 348        DrawEvent(currentEvent, fCenterTime, color);
 349
 350
 351
 352
 353
 354    }
 355    leg->Draw("same");
 356
 357 }
 358 //________________________________________________________________
 359
 360 void AliToyMCDrawer::DrawEvent(AliToyMCEvent *currentEvent, Double_t centerTime, Int_t color){
 361   if(!fDispHist) DrawGeometry();
 362
 363   for(Int_t iTrack = 0; iTrack < currentEvent->GetNumberOfTracks(); iTrack++){
 364
 365     //Double_t alpha = TMath::ATan2(currentEvent->GetTrack(iTrack)->GetSpacePoint(0)->GetY(),currentEvent->GetTrack(iTrack)->GetSpacePoint(0)->GetX());
 366     // if( (     (  !(abs(alpha)>phiCut && abs(alpha) < TMath::Pi() - phiCut)  ) )      )continue;
 367     // if(alpha<0) {
 368     //  continue;
 369     //  delete tempTrack;
 370     // }
 371
 372     //if(currentEvent->GetTrack(iTrack)->GetSpacePoint(0)->GetY()<0) alpha +=TMath::Pi();
 373     //if( (       (alpha>phiCut && alpha < TMath::Pi() - phiCut) || ( alpha>TMath::Pi() + phiCut && alpha < TMath::TwoPi() - phiCut))   )continue;
 374
 375     // if(alpha > 0) continue;
 376
 377     //std::cout << TMath::ASin(tempTrack->GetParameter()[2]) << std::endl;
 378     // Double_t x = currentEvent->GetX();
 379     // Double_t y = currentEvent->GetY();
 380       // Double_t hlxpar[6];
 381       // currentEvent->GetTrack(iTrack)->GetHelixParameters(hlxpar,0);
 382       // Double_t trackPhi = hlxpar[2];
 383       // Double_t phiCut = 1.47;
 384       // std::cout << "bef phi cut "<< trackPhi<< std::endl;
 385
 386       // std::cout << "aft phi cut " << std::endl;
 387       // Double_t trackEta = currentEvent->GetTrack(iTrack)->GetEta();
 388       //  Double_t z = currentEvent->GetZ();
 389
 390
 391
 392
 393
 394       const AliToyMCTrack *tempTrack = currentEvent->GetTrack(iTrack);
 395       //AliToyMCTrack *tempTrack = new AliToyMCTrack(*currentEvent->GetTrack(iTrack));
 396       DrawTrack(tempTrack, centerTime, currentEvent->GetT0(),color);
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407     }
 408
 409
 410
 411
 412 }
 413
 414 //________________________________________________________________
 415 void AliToyMCDrawer::DrawLaserEvent(Int_t nLaserEvents, Int_t side, Int_t rod, Int_t bundle, Int_t beam)
 416 {
 417   //
 418   //
 419   //
 420
 421   fPoints->Clear();
 422   fDistPoints->Clear();
 423   if (!ConnectInputTree()) return;
 424
 425   AliTPCLaserTrack::LoadTracks();
 426   TObjArray *arr=AliTPCLaserTrack::GetTracks();
 427
 428   DrawGeometry();
 429
 430   Int_t laserEvents=0;
 431   for (Int_t iev=0;iev<fInputTree->GetEntries();++iev) {
 432     fInputTree->GetEntry(iev);
 433     if (fEvent->GetEventType()!=AliToyMCEvent::kLaser) continue;
 434     if (fEvent->GetNumberOfTracks()!=arr->GetEntriesFast()) {
 435       AliError(Form("Wrong number of tracks in the laser event: %d!=%d",fEvent->GetNumberOfTracks(),arr->GetEntriesFast()));
 436       continue;
 437     }
 438
 439     for (Int_t iTrack=0; iTrack<fEvent->GetNumberOfTracks();++iTrack){
 440       AliTPCLaserTrack *laserTrack = (AliTPCLaserTrack*)arr->UncheckedAt(iTrack);
 441       if (side  >-1 && laserTrack->GetSide()   != side   ) continue;
 442       if (rod   >-1 && laserTrack->GetRod()    != rod    ) continue;
 443       if (bundle>-1 && laserTrack->GetBundle() != bundle ) continue;
 444       if (beam  >-1 && laserTrack->GetBeam()   != beam   ) continue;
 445       const AliToyMCTrack *track = fEvent->GetTrack(iTrack);
 446       DrawTrack(track,0,fEvent->GetT0(),kRed);
 447     }
 448
 449     ++laserEvents;
 450     if (laserEvents==nLaserEvents) break;
 451   }
 452
 453
 454 }
 455
 456 //________________________________________________________________
 457 void AliToyMCDrawer::DrawTrack(const AliToyMCTrack *track,  Double_t centerTime, Double_t currentEventTime, Int_t color){
 458   if(!fDispHist) DrawGeometry();
 459
 460   TPolyMarker3D *disttrackpoints = new((*fDistPoints)[fDistPoints->GetEntriesFast()]) TPolyMarker3D();
 461   TPolyMarker3D *trackpoints = new((*fPoints)[fPoints->GetEntriesFast()]) TPolyMarker3D();
 462
 463   Double_t currentEventTimeRelToCentral = centerTime - currentEventTime;
 464   Int_t nDistPoints = track->GetNumberOfDistSpacePoints();
 465   Int_t nPoints = track->GetNumberOfSpacePoints();
 466
 467   for(Int_t iITSPoint = 0; iITSPoint< track->GetNumberOfITSPoints();iITSPoint++){
 468
 469     Double_t xp = track->GetITSPoint(iITSPoint)->GetX();
 470     Double_t yp = track->GetITSPoint(iITSPoint)->GetY();
 471     Double_t zp = track->GetITSPoint(iITSPoint)->GetZ();
 472     Double_t zDrifted =  zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
 473     trackpoints->SetNextPoint(zDrifted,xp,yp);
 474   }
 475
 476   for(Int_t iPoint = 0; iPoint< (nDistPoints<nPoints?nPoints:nDistPoints);iPoint++){
 477     if(iPoint<nPoints) {
 478
 479
 480       Double_t xp = track->GetSpacePoint(iPoint)->GetX();
 481       Double_t yp = track->GetSpacePoint(iPoint)->GetY();
 482       Double_t zp = track->GetSpacePoint(iPoint)->GetZ();
 483       Double_t zDrifted =  zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
 484       //Double_t zDrifted = (zp/TMath::Abs(zp))*fMaxZ0  -(zp/TMath::Abs(zp))* fDriftVel*(track->GetSpacePoint(iPoint)->GetTimeBin()      -centerTime   );
 485       Float_t xyzp[3] = {xp,yp,zp};
 486       AliTrackPoint p;
 487       p.SetXYZ(xyzp);
 488       Float_t tempcov[6] = {0};
 489       p.SetCov(tempcov);
 490       Int_t sec = track->GetSpacePoint(iPoint)->GetDetector();
 491       Double_t angle=((sec%18)*20.+10.)/TMath::RadToDeg();
 492       AliTrackPoint prot = p.Rotate(-angle);
 493       xp = prot.GetX();
 494       yp = prot.GetY();
 495
 496       Double_t ztime=zDrifted;
 497       Double_t globx=xp;
 498       Double_t globy=yp;
 499
 500 //       if (fProje)
 501
 502
 503       if(track->GetSpacePoint(iPoint)->GetRow()!=255) {
 504         if(TMath::Abs(zDrifted)<fMaxZ0 && zDrifted*zp >=0 /*&&TMath::Sqrt(xp*xp + yp*yp)<fIFCRadius*/) trackpoints->SetNextPoint(zDrifted,xp,yp);
 505
 506       }
 507       else std::cout << "row == " << track->GetSpacePoint(iPoint)->GetRow() << std::endl;
 508
 509     }
 510     if(iPoint<nDistPoints) {
 511       Double_t xpdist = track->GetDistortedSpacePoint(iPoint)->GetX();
 512       Double_t ypdist = track->GetDistortedSpacePoint(iPoint)->GetY();
 513       Double_t zpdist = track->GetDistortedSpacePoint(iPoint)->GetZ();
 514       //std::cout << zpdist << std::endl;
 515
 516       Float_t xyzpdist[3] = {xpdist,ypdist,zpdist};
 517       AliTrackPoint pdist;
 518       pdist.SetXYZ(xyzpdist);
 519       Float_t tempcovdist[6] = {0};
 520       pdist.SetCov(tempcovdist);
 521       Int_t secdist = track->GetDistortedSpacePoint(iPoint)->GetDetector();
 522       Double_t angledist=((secdist%18)*20.+10.)/TMath::RadToDeg();
 523       AliTrackPoint protdist = pdist.Rotate(-angledist);
 524       xpdist = protdist.GetX();
 525       ypdist = protdist.GetY();
 526
 527
 528       UInt_t sector = track->GetDistortedSpacePoint(iPoint)->GetDetector();
 529       UInt_t row = track->GetDistortedSpacePoint(iPoint)->GetRow();
 530       UInt_t pad  = track->GetDistortedSpacePoint(iPoint)->GetPad();
 531
 532       Int_t nPads = fTPCParam->GetNPads(sector,row);
 533       Int_t intPad = TMath::Nint(Float_t(pad+nPads/2));
 534
 535       Float_t xyz[3];
 536       //std::cout <<"sector: " << sector << " row: " << row << " pad: " <<pad<< std::endl;
 537       fRoc->GetPositionGlobal(sector,row,intPad,xyz);
 538
 539       Double_t zDrifteddist = (zpdist/TMath::Abs(zpdist))*fMaxZ0  -(zpdist/TMath::Abs(zpdist))* fDriftVel*(track->GetDistortedSpacePoint(iPoint)->GetTimeBin()      -centerTime   );
 540       if(row!=255){
 541
 542         if(TMath::Abs(zDrifteddist)<fMaxZ0 && zDrifteddist*zpdist>=0 /*&&TMath::Sqrt(xpdist*xpdist + ypdist*ypdist)<fIFCRadius*/) {
 543
 544           disttrackpoints->SetNextPoint(zDrifteddist,xpdist,ypdist);
 545           //if(TMath::Sqrt(xpdist*xpdist + ypdist*ypdist)<fIFCRadius) std::cout << "fMaxZ0 " << fMaxZ0 <<" inside " << xpdist << " " << "zpdist "  << zpdist << " " << "zDrifteddist "<< zDrifteddist << " " << zDrifteddist*zpdist << std::endl;
 546         }
 547       }
 548       else std::cout << "row == " << row << std::endl;
 549       //Double_t zDrifteddist =(zpdist/TMath::Abs(zpdist))*fMaxZ0  -(zpdist/TMath::Abs(zpdist))*(currentEvent->GetTrack(iTrack)->GetDistortedSpacePoint(iPoint)->GetTimeBin()- currentEvent->GetT0() )* fDriftVel;
 550
 551
 552     }
 553     //  if( (       (trackPhi>phiCut && trackPhi < TMath::Pi() - phiCut) || ( trackPhi>TMath::Pi() + phiCut && trackPhi < TMath::TwoPi() - phiCut))   ) {
 554
 555
 556
 557
 558
 559   }
 560   if(0){
 561     for(Int_t iTRDPoint = 0; iTRDPoint< track->GetNumberOfTRDPoints();iTRDPoint++){
 562
 563       Double_t xp = track->GetTRDPoint(iTRDPoint)->GetX();
 564       Double_t yp = track->GetTRDPoint(iTRDPoint)->GetY();
 565       Double_t zp = track->GetTRDPoint(iTRDPoint)->GetZ();
 566       Double_t zDrifted =  zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
 567       trackpoints->SetNextPoint(zDrifted,xp,yp);
 568   }
 569 }
 570   if(1){
 571     if(trackpoints && trackpoints->GetN()>0) {
 572       //   trackpoints->SetMarkerColor(1+currentEvent->GetEventNumber()%9);
 573       //trackpoints->SetMarkerStyle(6);
 574       trackpoints->Draw("same");
 575     }
 576     if(disttrackpoints && disttrackpoints->GetN()>0) {
 577       //
 578       //disttrackpoints->SetMarkerStyle(6);
 579       disttrackpoints->SetMarkerColor(color);
 580
 581
 582       disttrackpoints->Draw("same");
 583
 584     }
 585
 586       }
 587
 588
 589
 590
 591 }
 592
 593 //________________________________________________________________
 594 void AliToyMCDrawer::DrawGeometry() {
 595
 596
 597   //delete fDispGraph;
 598   //fDispGraph = new TGraph2D();
 599   delete fDispHist;
 600
 601   Double_t timeZmin=-fMaxZ0;
 602   Double_t timeZmax= fMaxZ0;
 603   Double_t globXmin=-(fOFCRadius +10);
 604   Double_t globXmax=  fOFCRadius +10 ;
 605   Double_t globYmin=-(fOFCRadius +10);
 606   Double_t globYmax=  fOFCRadius +10 ;
 607
 608   const Double_t epsilon=.001;
 609
 610   TString title;
 611   if (fTimeZmax>fTimeZmin) {
 612     timeZmin=fTimeZmin;
 613     timeZmax=fTimeZmax;
 614   }
 615   if (fGlobalXmax>fGlobalXmin) {
 616     globXmin=fGlobalXmin;
 617     globXmax=fGlobalXmax;
 618   }
 619   if (fGlobalYmax>fGlobalYmin) {
 620     globYmin=fGlobalYmin;
 621     globYmax=fGlobalYmax;
 622   }
 623   if (fProjectionType=="XYT"){
 624     fDispHist = new TH3F("fDispHist",";#it{z} [cm]; #it{x} [cm]; #it{y} [cm]",
 625                          100, timeZmin-10, timeZmax+10,
 626                          100, globXmin, globXmax ,
 627                          100, globYmin, globYmax);
 628   } else if (fProjectionType=="XT"){
 629     fDispHist = new TH2F("fDispHist",";#it{z} [cm]; #it{x} [cm]",
 630                          100, timeZmin-10, timeZmax+10,
 631                          100, globXmin, globXmax);
 632   } else if (fProjectionType=="YT"){
 633     fDispHist = new TH2F("fDispHist",";#it{z} [cm]; #it{y} [cm]",
 634                          100, timeZmin-10, timeZmax+10,
 635                          100, globYmin, globYmax);
 636   } else if (fProjectionType=="RT"){
 637     fDispHist = new TH2F("fDispHist",";#it{z} [cm]; #it{r} [cm]",
 638                          100, timeZmin-10, timeZmax+10,
 639                          100, globYmin, globYmax);
 640   } else {
 641     AliError(Form("Display Format not known: %s",fProjectionType.Data()));
 642     return;
 643   }
 644
 645
 646
 647   //if(!fDispGraph) fDispGraph = new TGraph();
 648
 649   //fDispGraph->Clear();
 650   fDispHist->SetStats(0);
 651   fDispHist->Draw();
 652   gPad->SetPhi(0);
 653   gPad->SetTheta(0);
 654
 655   TPolyLine3D *endCap1 = 0x0;
 656   if (timeZmin-epsilon<-fMaxZ0)
 657     endCap1=new TPolyLine3D();
 658   printf("time: %p, %.2f,, %.2f, %.2f, %.2f, %d\n",endCap1,fTimeZmin,fTimeZmax,timeZmin-epsilon,fMaxZ0, timeZmin-epsilon<-fMaxZ0);
 659   TPolyLine3D *endCap2 = 0x0;
 660   if (timeZmax+epsilon> fMaxZ0)
 661     endCap2=new TPolyLine3D();
 662
 663   TPolyLine3D *outerCE = 0x0;
 664   if (timeZmin<0 && timeZmax>0 )
 665     outerCE=new TPolyLine3D();
 666
 667   TPolyLine3D *cage[18] ={0x0};
 668
 669   TPolyLine3D *innerCage[18] ={0x0};
 670
 671   TPolyLine3D *innerEndCap1 = 0x0;
 672   if (timeZmin-epsilon<-fMaxZ0)
 673     innerEndCap1=new TPolyLine3D();
 674
 675   TPolyLine3D *innerEndCap2 = 0x0;
 676   if (timeZmax+epsilon> fMaxZ0)
 677     innerEndCap2=new TPolyLine3D();
 678
 679   TPolyLine3D *innerCE = 0x0;
 680   if (timeZmin<0 && timeZmax>0 )
 681     innerCE=new TPolyLine3D();
 682
 683   Int_t iPoint=0;
 684   Double_t angle    = 0.;
 685   Double_t globalX  = 0.;
 686   Double_t globalY  = 0.;
 687   Double_t globalXi = 0.;
 688   Double_t globalYi = 0.;
 689
 690   for(Int_t i = 0; i<18; i++){
 691     angle    = i*TMath::TwoPi()/18;
 692     globalX  = fOFCRadius*TMath::Cos(angle);
 693     globalY  = fOFCRadius*TMath::Sin(angle);
 694     globalXi = fIFCRadius*TMath::Cos(angle);
 695     globalYi = fIFCRadius*TMath::Sin(angle);
 696
 697     cage[iPoint] = new TPolyLine3D();
 698     cage[iPoint]->SetPoint(0,timeZmin,globalX ,globalY) ;
 699     cage[iPoint]->SetPoint(1, timeZmax,globalX ,globalY) ;
 700     innerCage[iPoint] = new TPolyLine3D();
 701     innerCage[iPoint]->SetPoint(0,timeZmin,globalXi ,globalYi) ;
 702     innerCage[iPoint]->SetPoint(1, timeZmax,globalXi ,globalYi) ;
 703
 704     // only draw if inside range
 705     if (endCap1) { endCap1->SetPoint(i,timeZmax, globalX, globalY); }
 706     if (endCap2) { endCap2->SetPoint(i,timeZmin, globalX, globalY); }
 707     if (outerCE) { outerCE->SetPoint(i,      0., globalX, globalY); }
 708
 709     if (innerEndCap1) { innerEndCap1->SetPoint(i, timeZmax, globalXi, globalYi); }
 710     if (innerEndCap2) { innerEndCap2->SetPoint(i, timeZmin, globalXi, globalYi); }
 711     if (innerCE)      {      innerCE->SetPoint(i,       0., globalXi, globalYi); }
 712
 713     innerCage[iPoint]->Draw("same");
 714
 715     if(!(i%2))
 716       cage[iPoint]->Draw("same");
 717
 718     ++iPoint;
 719   }
 720
 721   //
 722   // close endplate and CE polygons
 723   //
 724   Int_t i=18;
 725   angle    = i*TMath::TwoPi()/18;
 726   globalX  = fOFCRadius*TMath::Cos(angle);
 727   globalY  = fOFCRadius*TMath::Sin(angle);
 728   globalXi = fIFCRadius*TMath::Cos(angle);
 729   globalYi = fIFCRadius*TMath::Sin(angle);
 730
 731   // only draw if inside range
 732   if (endCap1) { endCap1->SetPoint(i, timeZmax, globalX, globalY); }
 733   if (endCap2) { endCap2->SetPoint(i, timeZmin, globalX, globalY); }
 734   if (outerCE) { outerCE->SetPoint(i,       0., globalX, globalY); }
 735
 736   if (innerEndCap1) { innerEndCap1->SetPoint(i, timeZmax, globalXi, globalYi); }
 737   if (innerEndCap2) { innerEndCap2->SetPoint(i, timeZmin, globalXi, globalYi); }
 738   if (innerCE)      { innerCE     ->SetPoint(i,       0., globalXi, globalYi); }
 739
 740   if (endCap1) { endCap1->Draw("same"); }
 741   if (endCap2) { endCap2->Draw("same"); }
 742   if (outerCE) { outerCE->Draw("same"); }
 743
 744   if (innerEndCap1) { innerEndCap1->Draw("same"); }
 745   if (innerEndCap2) { innerEndCap2->Draw("same"); }
 746   if (innerCE)      { innerCE     ->Draw("same");      }
 747
 748 }
 749
 750 //________________________________________________________________
 751 Bool_t AliToyMCDrawer::ConnectInputTree()
 752 {
 753   //
 754   //
 755   //
 756
 757   if(fFileName.IsNull()) {
 758     AliError("no input file provided, using default (toyMC.root)");
 759     fFileName = "toyMC.root";
 760   }
 761   if (fInFile && fInFile->GetName()==fFileName) return kTRUE;
 762   delete fInFile;
 763   fInFile=0x0;
 764   delete fInputTree;
 765   fInputTree=0x0;
 766
 767   fInFile = new TFile(fFileName.Data(),"read");
 768   if (!fInFile || !fInFile->IsOpen() || fInFile->IsZombie() ) {
 769     delete fInFile;
 770     return kFALSE;
 771   }
 772
 773   gROOT->cd();
 774
 775   fInputTree = dynamic_cast<TTree*> (fInFile->Get("toyMCtree"));
 776
 777   if (!fInputTree){
 778     AliError("Could not connect tree!\n");
 779     delete fInFile;
 780     fInFile=0x0;
 781     return kFALSE;
 782   }
 783
 784   fInputTree->SetBranchAddress("event",&fEvent);
 785   return kTRUE;
 786 }