10 #include <TClonesArray.h>
11 #include <TPolyLine3D.h>
12 #include <TPolyLine.h>
13 #include <TPolyMarker3D.h>
14 #include <TPolyMarker.h>
18 #include <TGeoGlobalMagField.h>
22 #include <AliCDBManager.h>
23 #include <AliTPCParam.h>
24 #include <AliGeomManager.h>
25 #include <AliTPCcalibDB.h>
26 #include <AliTrackPointArray.h>
27 #include <AliCluster.h>
29 #include <AliTPCLaserTrack.h>
32 #include "AliToyMCDrawer.h"
33 #include "AliToyMCEvent.h"
34 #include "AliToyMCTrack.h"
36 // Visualization class. To use
38 // AliToyMCDrawer* draw = new AliToyMCDrawer()
39 // draw->SetFileName("path/to/toyMC.root")
41 // draw->FillEventArray(Int_t centerEventNumber)
43 // draw->FillEventArray(Double_t time)
44 // to display with a certain event in the center or at a certain time
46 // draw->DrawEvents(Bool_t both, Bool_t before)
47 // where "both" will display events before and after the middle event and
48 // before will show also events before (after) the middle event if true (false)
49 // when "both" is false
52 ClassImp(AliToyMCDrawer);
54 AliToyMCDrawer::AliToyMCDrawer()
69 ,fRoc(AliTPCROC::Instance())
72 ,fProjectionType("XYT")
80 fEventArray = new TClonesArray("AliToyMCEvent");
82 fTPCParam = AliTPCcalibDB::Instance()->GetParameters();
83 fTPCParam->ReadGeoMatrices();
84 fDriftVel = fTPCParam->GetDriftV();
85 fMaxZ0 =fTPCParam->GetZLength();
86 fTimeRange = 2*fMaxZ0/fDriftVel;
88 //________________________________________________________________
89 AliToyMCDrawer::AliToyMCDrawer(const AliToyMCDrawer &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)
105 ,fPoints(drawer.fPoints)
106 ,fDistPoints(drawer.fDistPoints)
107 ,fProjectionType("XYT")
117 //_____________________________________________________
118 AliToyMCDrawer& AliToyMCDrawer::operator = (const AliToyMCDrawer &drawer)
120 //assignment operator
121 if (&drawer == this) return *this;
122 new (this) AliToyMCDrawer(drawer);
127 //________________________________________________________________
128 AliToyMCDrawer::~AliToyMCDrawer()
137 //________________________________________________________________
138 Int_t AliToyMCDrawer::FillEventArray(Double_t snapShotTime)
140 if(fFileName.IsNull()) {
141 std::cout << "no input file provided, using default (toyMC.root)" << std::endl;
142 fFileName = "toyMC.root";
145 fInFile = new TFile(fFileName.Data(),"read");
147 fInputTree = dynamic_cast<TTree*> (fInFile->Get("toyMCtree"));
148 fInputTree->SetBranchAddress("event",&fEvent);
149 fEventArray->Clear();
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;
158 if(leftSearchIndex<0) leftSearchIndex = 0;
160 //fCenterTime = snapShotTime;
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)
168 // printf("close from start, lefttime = %2.2f\n",leftTime);
169 firstEventIndex = leftSearchIndex+1;
170 while( TMath::Abs(leftTime- snapShotTime) < 5e-9){
172 // printf("close from start, lefttime = %2.2f\n",leftTime);
173 if(fInputTree->GetEvent(firstEventIndex-1))
175 leftTime=fEvent->GetT0();
176 // printf("close from start, lefttime = %2.2f\n",leftTime);
185 Int_t rightSearchIndex = leftSearchIndex;
188 Int_t direction = 1; //go right
189 if(leftTime > snapShotTime) {
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;
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;
209 rightTime = fEvent->GetT0();
210 printf("searching, leftime %f, righttim %f\n",leftTime,rightTime);
212 if (direction==-1) rightSearchIndex = leftSearchIndex;
213 firstEventIndex = rightSearchIndex;
216 // fInputTree->GetEvent(firstEventIndex);
217 //std::cout <<"first event after sn time: " << fEvent->GetT0() << std::endl;
218 return FillEventArray(firstEventIndex, snapShotTime);
221 //________________________________________________________________
222 Int_t AliToyMCDrawer::FillEventArray(Int_t middleEventNbr, Double_t snapShotTime)
224 if(fFileName.IsNull()) {
225 std::cout << "no input file provided, using default (toyMC.root)" << std::endl;
226 fFileName = "toyMC.root";
228 if(!fInFile) fInFile = new TFile(fFileName.Data(),"read");
232 fInputTree = dynamic_cast<TTree*> (fInFile->Get("toyMCtree"));
233 // fInputTree->Add(fileName);
234 fInputTree->SetBranchAddress("event",&fEvent);
236 Double_t centerEventTime;
237 if(fInputTree->GetEvent(middleEventNbr))
239 centerEventTime = fEvent->GetT0();
244 if(snapShotTime<0) fCenterTime = centerEventTime;
245 else fCenterTime = snapShotTime;
247 fEventArray->Clear();
248 if(fInputTree->GetEvent(middleEventNbr)){
249 new((*fEventArray)[fEventArray->GetEntriesFast()]) AliToyMCEvent(*fEvent);
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);
259 if(!fInputTree->GetEvent(eventIndex)) break;
260 currentTime = fEvent->GetT0();
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) {
269 new((*fEventArray)[fEventArray->GetEntriesFast()]) AliToyMCEvent(*fEvent);
271 if(!fInputTree->GetEvent(eventIndex)) break;
272 currentTime = fEvent->GetT0();
276 std::cout << " end fill" << std::endl;
281 printf("Selected event number (%d) out of range!\n",middleEventNbr);
286 //________________________________________________________________
287 void AliToyMCDrawer::DrawEvents(Bool_t both, Bool_t before)
292 fDistPoints->Clear();
297 Double_t phiCut = 1.;
298 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));
299 leg->AddEntry((TObject*)0,Form("%2.2f<#phi<#pi-%2.2f && #pi+%2.2f<#phi<2#pi-%2.2f",phiCut,phiCut,phiCut,phiCut),"");
304 for(Int_t iEvent = 0; iEvent<fEventArray->GetEntriesFast(); iEvent++){
307 AliToyMCEvent *currentEvent = static_cast<AliToyMCEvent*> (fEventArray->At(iEvent));
308 Double_t currentEventTime = currentEvent->GetT0();
310 if((1000000*currentEventTime)>((1000000*fCenterTime)+0.05 )) {
311 printf("larger iEvent: %d, current %f, center %f, ntracks: %d\n",iEvent, currentEventTime,fCenterTime, currentEvent->GetNumberOfTracks());
313 if ( !(both || !before)) continue;
316 if( currentEventTime<=fCenterTime)// && !(currentEventTime==fCenterTime &&both ))
318 printf("smaller iEvent: %d, current %f, center %f, ntracks: %d\n",iEvent, currentEventTime,fCenterTime, currentEvent->GetNumberOfTracks());
320 if ( !(both || before)) continue;
323 //std:: cout << "iev " << iEvent << " " << color << std::endl;
324 Double_t currentEventTimeRelToCentral = fCenterTime - currentEventTime;
326 Int_t color = ((currentEvent->GetEventNumber())%10);
327 if(color == 0||color ==1) {
328 color==0?color=2:color=3;
329 // if(iEvent==0) color = 8;
330 // if(iEvent==1) color = 9;
334 char bef[] = "before snapshot time";
335 char aft[] = "after snapshot time";
338 if( 1000000*(currentEventTime-fCenterTime) > 0.5) {when = aft;
342 TGraph *temp = new TGraph();
343 temp->SetName(Form("temp%d",currentEvent->GetEventNumber()));
344 temp->SetLineColor(color);
345 temp->SetLineWidth(10);
347 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");
350 DrawEvent(currentEvent, fCenterTime, color);
360 //________________________________________________________________
362 void AliToyMCDrawer::DrawEvent(AliToyMCEvent *currentEvent, Double_t centerTime, Int_t color){
363 if(!fDispHist) DrawGeometry();
365 for(Int_t iTrack = 0; iTrack < currentEvent->GetNumberOfTracks(); iTrack++){
367 //Double_t alpha = TMath::ATan2(currentEvent->GetTrack(iTrack)->GetSpacePoint(0)->GetY(),currentEvent->GetTrack(iTrack)->GetSpacePoint(0)->GetX());
368 // if( ( ( !(abs(alpha)>phiCut && abs(alpha) < TMath::Pi() - phiCut) ) ) )continue;
374 //if(currentEvent->GetTrack(iTrack)->GetSpacePoint(0)->GetY()<0) alpha +=TMath::Pi();
375 //if( ( (alpha>phiCut && alpha < TMath::Pi() - phiCut) || ( alpha>TMath::Pi() + phiCut && alpha < TMath::TwoPi() - phiCut)) )continue;
377 // if(alpha > 0) continue;
379 //std::cout << TMath::ASin(tempTrack->GetParameter()[2]) << std::endl;
380 // Double_t x = currentEvent->GetX();
381 // Double_t y = currentEvent->GetY();
382 // Double_t hlxpar[6];
383 // currentEvent->GetTrack(iTrack)->GetHelixParameters(hlxpar,0);
384 // Double_t trackPhi = hlxpar[2];
385 // Double_t phiCut = 1.47;
386 // std::cout << "bef phi cut "<< trackPhi<< std::endl;
388 // std::cout << "aft phi cut " << std::endl;
389 // Double_t trackEta = currentEvent->GetTrack(iTrack)->GetEta();
390 // Double_t z = currentEvent->GetZ();
396 const AliToyMCTrack *tempTrack = currentEvent->GetTrack(iTrack);
397 //AliToyMCTrack *tempTrack = new AliToyMCTrack(*currentEvent->GetTrack(iTrack));
398 if (fProjectionType.Length()==3){
399 DrawTrack(tempTrack, centerTime, currentEvent->GetT0(),color);
400 } else if (fProjectionType.Length()==2){
401 DrawTrack2D(tempTrack, centerTime, currentEvent->GetT0(),color);
407 //________________________________________________________________
408 void AliToyMCDrawer::DrawLaserEvent(Int_t nLaserEvents, Int_t side, Int_t rod, Int_t bundle, Int_t beam)
416 fDistPoints->Clear();
418 if (!ConnectInputTree()) return;
420 AliTPCLaserTrack::LoadTracks();
421 TObjArray *arr=AliTPCLaserTrack::GetTracks();
423 if (fProjectionType.Length()==3){
425 } else if (fProjectionType.Length()==2){
430 for (Int_t iev=0;iev<fInputTree->GetEntries();++iev) {
431 fInputTree->GetEntry(iev);
432 if (fEvent->GetEventType()!=AliToyMCEvent::kLaser) continue;
433 if (fEvent->GetNumberOfTracks()!=arr->GetEntriesFast()) {
434 AliError(Form("Wrong number of tracks in the laser event: %d!=%d",fEvent->GetNumberOfTracks(),arr->GetEntriesFast()));
438 for (Int_t iTrack=0; iTrack<fEvent->GetNumberOfTracks();++iTrack){
439 AliTPCLaserTrack *laserTrack = (AliTPCLaserTrack*)arr->UncheckedAt(iTrack);
440 if (side >-1 && laserTrack->GetSide() != side ) continue;
441 if (rod >-1 && laserTrack->GetRod() != rod ) continue;
442 if (bundle>-1 && laserTrack->GetBundle() != bundle ) continue;
443 if (beam >-1 && laserTrack->GetBeam() != beam ) continue;
444 const AliToyMCTrack *track = fEvent->GetTrack(iTrack);
445 if (fProjectionType.Length()==3){
446 DrawTrack(track,0,fEvent->GetT0(),kRed);
447 } else if (fProjectionType.Length()==2){
448 DrawTrack2D(track,0,fEvent->GetT0(),kRed);
454 if (laserEvents==nLaserEvents) break;
460 //________________________________________________________________
461 void AliToyMCDrawer::DrawTrack(const AliToyMCTrack *track, Double_t centerTime, Double_t currentEventTime, Int_t color){
462 if(!fDispHist) DrawGeometry();
464 fPoints = new TClonesArray("TPolyMarker3D");
465 fDistPoints = new TClonesArray("TPolyMarker3D");
468 TPolyMarker3D *disttrackpoints = new((*fDistPoints)[fDistPoints->GetEntriesFast()]) TPolyMarker3D();
469 TPolyMarker3D *trackpoints = new((*fPoints)[fPoints->GetEntriesFast()]) TPolyMarker3D();
471 Double_t currentEventTimeRelToCentral = centerTime - currentEventTime;
472 Int_t nDistPoints = track->GetNumberOfDistSpacePoints();
473 Int_t nPoints = track->GetNumberOfSpacePoints();
475 for(Int_t iITSPoint = 0; iITSPoint< track->GetNumberOfITSPoints();iITSPoint++){
477 Double_t xp = track->GetITSPoint(iITSPoint)->GetX();
478 Double_t yp = track->GetITSPoint(iITSPoint)->GetY();
479 Double_t zp = track->GetITSPoint(iITSPoint)->GetZ();
480 Double_t zDrifted = zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
481 trackpoints->SetNextPoint(zDrifted,xp,yp);
484 for(Int_t iPoint = 0; iPoint< (nDistPoints<nPoints?nPoints:nDistPoints);iPoint++){
488 Double_t xp = track->GetSpacePoint(iPoint)->GetX();
489 Double_t yp = track->GetSpacePoint(iPoint)->GetY();
490 Double_t zp = track->GetSpacePoint(iPoint)->GetZ();
491 Double_t zDrifted = zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
492 //Double_t zDrifted = (zp/TMath::Abs(zp))*fMaxZ0 -(zp/TMath::Abs(zp))* fDriftVel*(track->GetSpacePoint(iPoint)->GetTimeBin() -centerTime );
493 Float_t xyzp[3] = {xp,yp,zp};
496 Float_t tempcov[6] = {0};
498 Int_t sec = track->GetSpacePoint(iPoint)->GetDetector();
499 Double_t angle=((sec%18)*20.+10.)/TMath::RadToDeg();
500 AliTrackPoint prot = p.Rotate(-angle);
504 // Double_t ztime=zDrifted;
505 // Double_t globx=xp;
506 // Double_t globy=yp;
511 if(track->GetSpacePoint(iPoint)->GetRow()!=255) {
512 if(TMath::Abs(zDrifted)<fMaxZ0 && zDrifted*zp >=0 /*&&TMath::Sqrt(xp*xp + yp*yp)<fIFCRadius*/) trackpoints->SetNextPoint(zDrifted,xp,yp);
515 else std::cout << "row == " << track->GetSpacePoint(iPoint)->GetRow() << std::endl;
518 if(iPoint<nDistPoints) {
519 Double_t xpdist = track->GetDistortedSpacePoint(iPoint)->GetX();
520 Double_t ypdist = track->GetDistortedSpacePoint(iPoint)->GetY();
521 Double_t zpdist = track->GetDistortedSpacePoint(iPoint)->GetZ();
522 //std::cout << zpdist << std::endl;
524 Float_t xyzpdist[3] = {xpdist,ypdist,zpdist};
526 pdist.SetXYZ(xyzpdist);
527 Float_t tempcovdist[6] = {0};
528 pdist.SetCov(tempcovdist);
529 Int_t secdist = track->GetDistortedSpacePoint(iPoint)->GetDetector();
530 Double_t angledist=((secdist%18)*20.+10.)/TMath::RadToDeg();
531 AliTrackPoint protdist = pdist.Rotate(-angledist);
532 xpdist = protdist.GetX();
533 ypdist = protdist.GetY();
536 UInt_t sector = track->GetDistortedSpacePoint(iPoint)->GetDetector();
537 UInt_t row = track->GetDistortedSpacePoint(iPoint)->GetRow();
538 UInt_t pad = track->GetDistortedSpacePoint(iPoint)->GetPad();
540 Int_t nPads = fTPCParam->GetNPads(sector,row);
541 Int_t intPad = TMath::Nint(Float_t(pad+nPads/2));
544 //std::cout <<"sector: " << sector << " row: " << row << " pad: " <<pad<< std::endl;
545 fRoc->GetPositionGlobal(sector,row,intPad,xyz);
547 Double_t zDrifteddist = (zpdist/TMath::Abs(zpdist))*fMaxZ0 -(zpdist/TMath::Abs(zpdist))* fDriftVel*(track->GetDistortedSpacePoint(iPoint)->GetTimeBin() -centerTime );
550 if(TMath::Abs(zDrifteddist)<fMaxZ0 && zDrifteddist*zpdist>=0 /*&&TMath::Sqrt(xpdist*xpdist + ypdist*ypdist)<fIFCRadius*/) {
552 disttrackpoints->SetNextPoint(zDrifteddist,xpdist,ypdist);
553 //if(TMath::Sqrt(xpdist*xpdist + ypdist*ypdist)<fIFCRadius) std::cout << "fMaxZ0 " << fMaxZ0 <<" inside " << xpdist << " " << "zpdist " << zpdist << " " << "zDrifteddist "<< zDrifteddist << " " << zDrifteddist*zpdist << std::endl;
556 else std::cout << "row == " << row << std::endl;
557 //Double_t zDrifteddist =(zpdist/TMath::Abs(zpdist))*fMaxZ0 -(zpdist/TMath::Abs(zpdist))*(currentEvent->GetTrack(iTrack)->GetDistortedSpacePoint(iPoint)->GetTimeBin()- currentEvent->GetT0() )* fDriftVel;
561 // if( ( (trackPhi>phiCut && trackPhi < TMath::Pi() - phiCut) || ( trackPhi>TMath::Pi() + phiCut && trackPhi < TMath::TwoPi() - phiCut)) ) {
569 for(Int_t iTRDPoint = 0; iTRDPoint< track->GetNumberOfTRDPoints();iTRDPoint++){
571 Double_t xp = track->GetTRDPoint(iTRDPoint)->GetX();
572 Double_t yp = track->GetTRDPoint(iTRDPoint)->GetY();
573 Double_t zp = track->GetTRDPoint(iTRDPoint)->GetZ();
574 Double_t zDrifted = zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
575 trackpoints->SetNextPoint(zDrifted,xp,yp);
579 if(trackpoints && trackpoints->GetN()>0) {
580 // trackpoints->SetMarkerColor(1+currentEvent->GetEventNumber()%9);
581 //trackpoints->SetMarkerStyle(6);
582 trackpoints->Draw("same");
584 if(disttrackpoints && disttrackpoints->GetN()>0) {
586 //disttrackpoints->SetMarkerStyle(6);
587 disttrackpoints->SetMarkerColor(color);
590 disttrackpoints->Draw("same");
601 //________________________________________________________________
602 void AliToyMCDrawer::DrawTrack2D(const AliToyMCTrack *track, Double_t centerTime, Double_t currentEventTime, Int_t color){
603 if(!fDispHist) DrawGeometry2D();
605 fPoints = new TClonesArray("TPolyMarker");
606 fDistPoints = new TClonesArray("TPolyMarker");
609 Double_t timeZmin=-fMaxZ0;
610 Double_t timeZmax= fMaxZ0;
611 Double_t globXmin=-(fOFCRadius +10);
612 Double_t globXmax= fOFCRadius +10 ;
613 Double_t globYmin=-(fOFCRadius +10);
614 Double_t globYmax= fOFCRadius +10 ;
616 // const Double_t epsilon=.001;
619 if (fTimeZmax>fTimeZmin) {
623 if (fGlobalXmax>fGlobalXmin) {
624 globXmin=fGlobalXmin;
625 globXmax=fGlobalXmax;
627 if (fGlobalYmax>fGlobalYmin) {
628 globYmin=fGlobalYmin;
629 globYmax=fGlobalYmax;
632 TPolyMarker *disttrackpoints = new((*fDistPoints)[fDistPoints->GetEntriesFast()]) TPolyMarker();
633 TPolyMarker *trackpoints = new((*fPoints)[fPoints->GetEntriesFast()]) TPolyMarker();
635 Double_t currentEventTimeRelToCentral = centerTime - currentEventTime;
636 Int_t nDistPoints = track->GetNumberOfDistSpacePoints();
637 Int_t nPoints = track->GetNumberOfSpacePoints();
639 for(Int_t iITSPoint = 0; iITSPoint< track->GetNumberOfITSPoints();iITSPoint++){
641 Double_t xp = track->GetITSPoint(iITSPoint)->GetX();
642 Double_t yp = track->GetITSPoint(iITSPoint)->GetY();
643 Double_t zp = track->GetITSPoint(iITSPoint)->GetZ();
644 Double_t zDrifted = zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
646 if ( xp<globXmin || xp>globXmax || yp<globYmin || yp>globYmax || zDrifted<timeZmin || zDrifted>timeZmax ) continue;
650 if (fProjectionType=="XT"){
653 } else if (fProjectionType=="YT"){
656 } else if (fProjectionType=="RT"){
658 y=TMath::Sqrt(xp*xp+yp*yp);
659 } else if (fProjectionType=="XY" || fProjectionType=="YX") {
663 trackpoints->SetNextPoint(x,y);
666 for(Int_t iPoint = 0; iPoint< (nDistPoints<nPoints?nPoints:nDistPoints);iPoint++){
670 Double_t xp = track->GetSpacePoint(iPoint)->GetX();
671 Double_t yp = track->GetSpacePoint(iPoint)->GetY();
672 Double_t zp = track->GetSpacePoint(iPoint)->GetZ();
673 Double_t zDrifted = zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
675 //Double_t zDrifted = (zp/TMath::Abs(zp))*fMaxZ0 -(zp/TMath::Abs(zp))* fDriftVel*(track->GetSpacePoint(iPoint)->GetTimeBin() -centerTime );
676 Float_t xyzp[3] = {xp,yp,zp};
679 Float_t tempcov[6] = {0};
681 Int_t sec = track->GetSpacePoint(iPoint)->GetDetector();
682 Double_t angle=((sec%18)*20.+10.)/TMath::RadToDeg();
683 AliTrackPoint prot = p.Rotate(-angle);
687 if ( xp<globXmin || xp>globXmax || yp<globYmin || yp>globYmax || zDrifted<timeZmin || zDrifted>timeZmax ) continue;
689 // Double_t ztime=zDrifted;
690 // Double_t globx=xp;
691 // Double_t globy=yp;
695 if (fProjectionType=="XT"){
698 } else if (fProjectionType=="YT"){
701 } else if (fProjectionType=="RT"){
703 y=TMath::Sqrt(xp*xp+yp*yp);
704 } else if (fProjectionType=="XY" || fProjectionType=="YX") {
709 if(track->GetSpacePoint(iPoint)->GetRow()!=255) {
710 if(TMath::Abs(zDrifted)<fMaxZ0 && zDrifted*zp >=0 /*&&TMath::Sqrt(xp*xp + yp*yp)<fIFCRadius*/) trackpoints->SetNextPoint(x,y);
713 else std::cout << "row == " << track->GetSpacePoint(iPoint)->GetRow() << std::endl;
716 if(iPoint<nDistPoints) {
717 Double_t xpdist = track->GetDistortedSpacePoint(iPoint)->GetX();
718 Double_t ypdist = track->GetDistortedSpacePoint(iPoint)->GetY();
719 Double_t zpdist = track->GetDistortedSpacePoint(iPoint)->GetZ();
721 //std::cout << zpdist << std::endl;
723 Float_t xyzpdist[3] = {xpdist,ypdist,zpdist};
725 pdist.SetXYZ(xyzpdist);
726 Float_t tempcovdist[6] = {0};
727 pdist.SetCov(tempcovdist);
728 Int_t secdist = track->GetDistortedSpacePoint(iPoint)->GetDetector();
729 Double_t angledist=((secdist%18)*20.+10.)/TMath::RadToDeg();
730 AliTrackPoint protdist = pdist.Rotate(-angledist);
731 xpdist = protdist.GetX();
732 ypdist = protdist.GetY();
734 UInt_t sector = track->GetDistortedSpacePoint(iPoint)->GetDetector();
735 UInt_t row = track->GetDistortedSpacePoint(iPoint)->GetRow();
736 UInt_t pad = track->GetDistortedSpacePoint(iPoint)->GetPad();
738 Int_t nPads = fTPCParam->GetNPads(sector,row);
739 Int_t intPad = TMath::Nint(Float_t(pad+nPads/2));
742 //std::cout <<"sector: " << sector << " row: " << row << " pad: " <<pad<< std::endl;
743 fRoc->GetPositionGlobal(sector,row,intPad,xyz);
745 Double_t zDrifteddist = (zpdist/TMath::Abs(zpdist))*fMaxZ0 -(zpdist/TMath::Abs(zpdist))* fDriftVel*(track->GetDistortedSpacePoint(iPoint)->GetTimeBin() -centerTime );
747 if ( xpdist<globXmin || xpdist>globXmax || ypdist<globYmin || ypdist>globYmax || zDrifteddist<timeZmin || zDrifteddist>timeZmax ) continue;
751 if (fProjectionType=="XT"){
754 } else if (fProjectionType=="YT"){
757 } else if (fProjectionType=="RT"){
759 y=TMath::Sqrt(xpdist*xpdist+ypdist*ypdist);
760 } else if (fProjectionType=="XY" || fProjectionType=="YX") {
767 if(TMath::Abs(zDrifteddist)<fMaxZ0 && zDrifteddist*zpdist>=0 /*&&TMath::Sqrt(xpdist*xpdist + ypdist*ypdist)<fIFCRadius*/) {
769 disttrackpoints->SetNextPoint(x,y);
770 //if(TMath::Sqrt(xpdist*xpdist + ypdist*ypdist)<fIFCRadius) std::cout << "fMaxZ0 " << fMaxZ0 <<" inside " << xpdist << " " << "zpdist " << zpdist << " " << "zDrifteddist "<< zDrifteddist << " " << zDrifteddist*zpdist << std::endl;
773 else std::cout << "row == " << row << std::endl;
774 //Double_t zDrifteddist =(zpdist/TMath::Abs(zpdist))*fMaxZ0 -(zpdist/TMath::Abs(zpdist))*(currentEvent->GetTrack(iTrack)->GetDistortedSpacePoint(iPoint)->GetTimeBin()- currentEvent->GetT0() )* fDriftVel;
778 // if( ( (trackPhi>phiCut && trackPhi < TMath::Pi() - phiCut) || ( trackPhi>TMath::Pi() + phiCut && trackPhi < TMath::TwoPi() - phiCut)) ) {
786 for(Int_t iTRDPoint = 0; iTRDPoint< track->GetNumberOfTRDPoints();iTRDPoint++){
788 Double_t xp = track->GetTRDPoint(iTRDPoint)->GetX();
789 Double_t yp = track->GetTRDPoint(iTRDPoint)->GetY();
790 Double_t zp = track->GetTRDPoint(iTRDPoint)->GetZ();
791 Double_t zDrifted = zp+(zp/TMath::Abs(zp))*currentEventTimeRelToCentral * fDriftVel;
793 if ( xp<globXmin || xp>globXmax || yp<globYmin || yp>globYmax || zDrifted<timeZmin || zDrifted>timeZmax ) continue;
797 if (fProjectionType=="XT"){
800 } else if (fProjectionType=="YT"){
803 } else if (fProjectionType=="RT"){
805 y=TMath::Sqrt(xp*xp+yp*yp);
806 } else if (fProjectionType=="XY" || fProjectionType=="YX") {
810 trackpoints->SetNextPoint(x,y);
814 if(trackpoints && trackpoints->GetN()>0) {
815 // trackpoints->SetMarkerColor(1+currentEvent->GetEventNumber()%9);
816 //trackpoints->SetMarkerStyle(6);
817 trackpoints->Draw("same");
819 if(disttrackpoints && disttrackpoints->GetN()>0) {
821 //disttrackpoints->SetMarkerStyle(6);
822 disttrackpoints->SetMarkerColor(color);
825 disttrackpoints->Draw("same");
836 //________________________________________________________________
837 void AliToyMCDrawer::DrawGeometry() {
841 //fDispGraph = new TGraph2D();
844 Double_t timeZmin=-fMaxZ0;
845 Double_t timeZmax= fMaxZ0;
846 Double_t globXmin=-(fOFCRadius +10);
847 Double_t globXmax= fOFCRadius +10 ;
848 Double_t globYmin=-(fOFCRadius +10);
849 Double_t globYmax= fOFCRadius +10 ;
851 const Double_t epsilon=.001;
854 if (fTimeZmax>fTimeZmin) {
858 if (fGlobalXmax>fGlobalXmin) {
859 globXmin=fGlobalXmin;
860 globXmax=fGlobalXmax;
862 if (fGlobalYmax>fGlobalYmin) {
863 globYmin=fGlobalYmin;
864 globYmax=fGlobalYmax;
866 fDispHist = new TH3F("fDispHist",";#it{z} (cm); #it{x} (cm); #it{y} (cm)",
867 100, timeZmin-10, timeZmax+10,
868 100, globXmin, globXmax ,
869 100, globYmin, globYmax);
871 //if(!fDispGraph) fDispGraph = new TGraph();
873 //fDispGraph->Clear();
874 fDispHist->SetStats(0);
879 TPolyLine3D *endCap1 = 0x0;
880 if (timeZmin-epsilon<-fMaxZ0)
881 endCap1=new TPolyLine3D();
882 printf("time: %p, %.2f,, %.2f, %.2f, %.2f, %d\n",endCap1,fTimeZmin,fTimeZmax,timeZmin-epsilon,fMaxZ0, timeZmin-epsilon<-fMaxZ0);
883 TPolyLine3D *endCap2 = 0x0;
884 if (timeZmax+epsilon> fMaxZ0)
885 endCap2=new TPolyLine3D();
887 TPolyLine3D *outerCE = 0x0;
888 if (timeZmin<0 && timeZmax>0 )
889 outerCE=new TPolyLine3D();
891 TPolyLine3D *cage[18] ={0x0};
893 TPolyLine3D *innerCage[18] ={0x0};
895 TPolyLine3D *innerEndCap1 = 0x0;
896 if (timeZmin-epsilon<-fMaxZ0)
897 innerEndCap1=new TPolyLine3D();
899 TPolyLine3D *innerEndCap2 = 0x0;
900 if (timeZmax+epsilon> fMaxZ0)
901 innerEndCap2=new TPolyLine3D();
903 TPolyLine3D *innerCE = 0x0;
904 if (timeZmin<0 && timeZmax>0 )
905 innerCE=new TPolyLine3D();
909 Double_t globalX = 0.;
910 Double_t globalY = 0.;
911 Double_t globalXi = 0.;
912 Double_t globalYi = 0.;
914 for(Int_t i = 0; i<18; i++){
915 angle = i*TMath::TwoPi()/18;
916 globalX = fOFCRadius*TMath::Cos(angle);
917 globalY = fOFCRadius*TMath::Sin(angle);
918 globalXi = fIFCRadius*TMath::Cos(angle);
919 globalYi = fIFCRadius*TMath::Sin(angle);
921 cage[iPoint] = new TPolyLine3D();
922 cage[iPoint]->SetPoint(0,timeZmin,globalX ,globalY) ;
923 cage[iPoint]->SetPoint(1, timeZmax,globalX ,globalY) ;
924 innerCage[iPoint] = new TPolyLine3D();
925 innerCage[iPoint]->SetPoint(0,timeZmin,globalXi ,globalYi) ;
926 innerCage[iPoint]->SetPoint(1, timeZmax,globalXi ,globalYi) ;
928 // only draw if inside range
929 if (endCap1) { endCap1->SetPoint(i,timeZmax, globalX, globalY); }
930 if (endCap2) { endCap2->SetPoint(i,timeZmin, globalX, globalY); }
931 if (outerCE) { outerCE->SetPoint(i, 0., globalX, globalY); }
933 if (innerEndCap1) { innerEndCap1->SetPoint(i, timeZmax, globalXi, globalYi); }
934 if (innerEndCap2) { innerEndCap2->SetPoint(i, timeZmin, globalXi, globalYi); }
935 if (innerCE) { innerCE->SetPoint(i, 0., globalXi, globalYi); }
937 innerCage[iPoint]->Draw("same");
940 cage[iPoint]->Draw("same");
946 // close endplate and CE polygons
949 angle = i*TMath::TwoPi()/18;
950 globalX = fOFCRadius*TMath::Cos(angle);
951 globalY = fOFCRadius*TMath::Sin(angle);
952 globalXi = fIFCRadius*TMath::Cos(angle);
953 globalYi = fIFCRadius*TMath::Sin(angle);
955 // only draw if inside range
956 if (endCap1) { endCap1->SetPoint(i, timeZmax, globalX, globalY); }
957 if (endCap2) { endCap2->SetPoint(i, timeZmin, globalX, globalY); }
958 if (outerCE) { outerCE->SetPoint(i, 0., globalX, globalY); }
960 if (innerEndCap1) { innerEndCap1->SetPoint(i, timeZmax, globalXi, globalYi); }
961 if (innerEndCap2) { innerEndCap2->SetPoint(i, timeZmin, globalXi, globalYi); }
962 if (innerCE) { innerCE ->SetPoint(i, 0., globalXi, globalYi); }
964 if (endCap1) { endCap1->Draw("same"); }
965 if (endCap2) { endCap2->Draw("same"); }
966 if (outerCE) { outerCE->Draw("same"); }
968 if (innerEndCap1) { innerEndCap1->Draw("same"); }
969 if (innerEndCap2) { innerEndCap2->Draw("same"); }
970 if (innerCE) { innerCE ->Draw("same"); }
974 //________________________________________________________________
975 void AliToyMCDrawer::DrawGeometry2D() {
979 //fDispGraph = new TGraph2D();
982 Double_t timeZmin=-fMaxZ0;
983 Double_t timeZmax= fMaxZ0;
984 Double_t globXmin=-(fOFCRadius +10);
985 Double_t globXmax= fOFCRadius +10 ;
986 Double_t globYmin=-(fOFCRadius +10);
987 Double_t globYmax= fOFCRadius +10 ;
989 const Double_t epsilon=.001;
992 if (fTimeZmax>fTimeZmin) {
996 if (fGlobalXmax>fGlobalXmin) {
997 globXmin=fGlobalXmin;
998 globXmax=fGlobalXmax;
1000 if (fGlobalYmax>fGlobalYmin) {
1001 globYmin=fGlobalYmin;
1002 globYmax=fGlobalYmax;
1004 TCanvas *c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject("cDrawer");
1006 if (fProjectionType=="XT"){
1007 if (!c) c=new TCanvas("cDrawer","Toy Drawer");
1008 fDispHist = new TH2F("fDispHist",";#it{z} (cm); #it{x} (cm)",
1009 100, timeZmin-10, timeZmax+10,
1010 100, globXmin, globXmax);
1011 } else if (fProjectionType=="YT"){
1012 if (!c) c=new TCanvas("cDrawer","Toy Drawer");
1013 fDispHist = new TH2F("fDispHist",";#it{z} (cm); #it{y} (cm)",
1014 100, timeZmin-10, timeZmax+10,
1015 100, globYmin, globYmax);
1016 } else if (fProjectionType=="RT"){
1017 if (!c) c=new TCanvas("cDrawer","Toy Drawer");
1018 fDispHist = new TH2F("fDispHist",";#it{z} (cm); #it{r} (cm)",
1019 100, timeZmin-10, timeZmax+10,
1020 100, globYmin, globYmax);
1021 } else if (fProjectionType=="YX"||fProjectionType=="XY"){
1022 if (!c) c=new TCanvas("cDrawer","Toy Drawer",800,800);
1023 c->SetLeftMargin(0.15);
1024 c->SetRightMargin(0.05);
1025 fDispHist = new TH2F("fDispHist",";#it{x} (cm); #it{y} (cm)",
1026 100, globXmin, globXmax,
1027 100, globYmin, globYmax);
1028 fDispHist->GetYaxis()->SetTitleOffset(1.5);
1030 AliError(Form("Display Format not known: %s",fProjectionType.Data()));
1036 //if(!fDispGraph) fDispGraph = new TGraph();
1038 //fDispGraph->Clear();
1039 fDispHist->SetStats(0);
1044 TPolyLine *endCap1 = 0x0;
1045 if (timeZmin-epsilon<-fMaxZ0)
1046 endCap1=new TPolyLine();
1047 printf("time: %p, %.2f,, %.2f, %.2f, %.2f, %d\n",endCap1,fTimeZmin,fTimeZmax,timeZmin-epsilon,fMaxZ0, timeZmin-epsilon<-fMaxZ0);
1048 TPolyLine *endCap2 = 0x0;
1049 if (timeZmax+epsilon> fMaxZ0)
1050 endCap2=new TPolyLine();
1052 TPolyLine *outerCE = 0x0;
1053 if (timeZmin<0 && timeZmax>0 )
1054 outerCE=new TPolyLine();
1056 TPolyLine *cage[18] ={0x0};
1058 TPolyLine *innerCage[18] ={0x0};
1060 TPolyLine *innerEndCap1 = 0x0;
1061 if (timeZmin-epsilon<-fMaxZ0)
1062 innerEndCap1=new TPolyLine();
1064 TPolyLine *innerEndCap2 = 0x0;
1065 if (timeZmax+epsilon> fMaxZ0)
1066 innerEndCap2=new TPolyLine();
1068 TPolyLine *innerCE = 0x0;
1069 if (timeZmin<0 && timeZmax>0 )
1070 innerCE=new TPolyLine();
1073 Double_t angle = 0.;
1074 Double_t globalX = 0.;
1075 Double_t globalY = 0.;
1076 Double_t globalXi = 0.;
1077 Double_t globalYi = 0.;
1078 Double_t globalXi2= 0.;
1079 Double_t globalYi2= 0.;
1082 if (fProjectionType=="YX"||fProjectionType=="XY") {
1083 for(Int_t i = 0; i<18; i++){
1084 angle = i*TMath::TwoPi()/18;
1085 globalX = fOFCRadius*TMath::Cos(angle);
1086 globalY = fOFCRadius*TMath::Sin(angle);
1087 globalXi = fIFCRadius*TMath::Cos(angle);
1088 globalYi = fIFCRadius*TMath::Sin(angle);
1089 globalXi = fIFCRadius*TMath::Cos(angle);
1090 globalYi = fIFCRadius*TMath::Sin(angle);
1091 globalXi2= 135.5*TMath::Cos(angle);
1092 globalYi2= 135.5*TMath::Sin(angle);
1094 if ( globalX<globXmin || globalX>globXmax || globalY<globYmin || globalY>globYmax ) {
1097 // in xy, abuse for sector boundaries
1098 cage[iPoint] = new TPolyLine();
1099 cage[iPoint]->SetPoint(0, globalX ,globalY) ;
1100 cage[iPoint]->SetPoint(1, globalXi ,globalYi) ;
1101 // only draw if inside range
1102 if (endCap1) { endCap1->SetPoint(iPoint, globalX, globalY); }
1104 if (innerEndCap1) { innerEndCap1->SetPoint(iPoint, globalXi, globalYi); }
1105 if (innerEndCap2) { innerEndCap2->SetPoint(iPoint, globalXi2, globalYi2); }
1107 cage[iPoint]->Draw("same");
1113 angle = 5*TMath::TwoPi()/18;
1114 globalX = fOFCRadius*TMath::Cos(angle);
1115 globalY = fOFCRadius*TMath::Sin(angle);
1116 globalXi = fIFCRadius*TMath::Cos(angle);
1117 globalYi = fIFCRadius*TMath::Sin(angle);
1118 globalXi = fIFCRadius*TMath::Cos(angle);
1119 globalYi = fIFCRadius*TMath::Sin(angle);
1120 globalXi2= 135.5*TMath::Cos(angle);
1121 globalYi2= 135.5*TMath::Sin(angle);
1122 if (endCap1) { endCap1->SetPoint(iPoint, 0, globalY); }
1124 if (innerEndCap1) { innerEndCap1->SetPoint(iPoint, 0, globalYi); }
1125 if (innerEndCap2) { innerEndCap2->SetPoint(iPoint, 0, globalYi2); }
1128 // close endplate and CE polygons
1131 angle = i*TMath::TwoPi()/18;
1132 globalX = fOFCRadius*TMath::Cos(angle);
1133 globalY = fOFCRadius*TMath::Sin(angle);
1134 globalXi = fIFCRadius*TMath::Cos(angle);
1135 globalYi = fIFCRadius*TMath::Sin(angle);
1136 globalXi2= 135.5*TMath::Cos(angle);
1137 globalYi2= 135.5*TMath::Sin(angle);
1139 // only draw if inside range
1140 if ( !(globalX<globXmin || globalX>globXmax || globalY<globYmin || globalY>globYmax) ) {
1141 if (endCap1) { endCap1->SetPoint(i, globalX, globalY); }
1142 // if (endCap2) { endCap2->SetPoint(i, globalX, globalY); }
1143 // if (outerCE) { outerCE->SetPoint(i, globalX, globalY); }
1145 if (innerEndCap1) { innerEndCap1->SetPoint(i, globalXi, globalYi); }
1146 if (innerEndCap2) { innerEndCap2->SetPoint(i, globalXi2, globalYi2); }
1147 // if (innerCE) { innerCE ->SetPoint(i, globalXi, globalYi); }
1151 if (endCap1) { endCap1->Draw("same"); }
1152 // if (endCap2) { endCap2->Draw("same"); }
1153 // if (outerCE) { outerCE->Draw("same"); }
1155 if (innerEndCap1) { innerEndCap1->Draw("same"); }
1156 if (innerEndCap2) { innerEndCap2->Draw("same"); }
1157 // if (innerCE) { innerCE ->Draw("same"); }
1160 //________________________________________________________________
1161 Bool_t AliToyMCDrawer::ConnectInputTree()
1167 if(fFileName.IsNull()) {
1168 AliError("no input file provided, using default (toyMC.root)");
1169 fFileName = "toyMC.root";
1171 if (fInFile && fInFile->GetName()==fFileName) return kTRUE;
1177 fInFile = new TFile(fFileName.Data(),"read");
1178 if (!fInFile || !fInFile->IsOpen() || fInFile->IsZombie() ) {
1185 fInputTree = dynamic_cast<TTree*> (fInFile->Get("toyMCtree"));
1188 AliError("Could not connect tree!\n");
1194 fInputTree->SetBranchAddress("event",&fEvent);