/************************************************************************** * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ ///////////////////////////////////////////////////////////// // // AliAnalysisTask to extract from ESD tracks the AliTrackPointArrays // with ITS points for selected tracks. This are the input data for alignment // // Author: A.Dainese, andrea.dainese@pd.infn.it ///////////////////////////////////////////////////////////// #include #include #include #include #include #include #include #include #include #include "AliLog.h" #include "AliGeomManager.h" #include "AliITSReconstructor.h" #include "AliITSgeomTGeo.h" #include "AliTrackPointArray.h" #include "AliESDInputHandler.h" #include "AliESDVertex.h" #include "AliESDtrack.h" #include "AliESDEvent.h" #include "AliESDfriend.h" #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" #include "AliAlignmentDataFilterITS.h" ClassImp(AliAlignmentDataFilterITS) //________________________________________________________________________ AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name): AliAnalysisTask(name,"task"), fOnlySPDFO(kFALSE), fGeometryFileName("geometry.root"), fITSRecoParam(0), fESD(0), fESDfriend(0), fListOfHistos(0), fspTree(0), fHistNevents(0), fHistNpoints(0), fHistPt(0), fHistLayer0(0), fHistLayer1(0), fHistLayer2(0), fHistLayer3(0), fHistLayer4(0), fHistLayer5(0), fntExtra(0), fntCosmicMatching(0) { // Constructor // Define input and output slots here // Input slot #0 works with a TChain DefineInput(0, TChain::Class()); // Output slot #0 writes into a TTree DefineOutput(0,TTree::Class()); //My private output // Output slot #1 writes into a TList DefineOutput(1,TList::Class()); //My private output } //________________________________________________________________________ AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS() { // Destructor if (fListOfHistos) { delete fListOfHistos; fListOfHistos = 0; } if (fspTree) { delete fspTree; fspTree = 0; } if (fHistNevents) { delete fHistNevents; fHistNevents = 0; } if (fHistNpoints) { delete fHistNpoints; fHistNpoints = 0; } if (fHistPt) { delete fHistPt; fHistPt = 0; } if (fHistLayer0) { delete fHistLayer0; fHistLayer0 = 0; } if (fHistLayer1) { delete fHistLayer1; fHistLayer1 = 0; } if (fHistLayer2) { delete fHistLayer2; fHistLayer2 = 0; } if (fHistLayer3) { delete fHistLayer3; fHistLayer3 = 0; } if (fHistLayer4) { delete fHistLayer4; fHistLayer4 = 0; } if (fHistLayer5) { delete fHistLayer5; fHistLayer5 = 0; } if (fntExtra) { delete fntExtra; fntExtra = 0; } if (fntCosmicMatching) { delete fntCosmicMatching; fntCosmicMatching = 0; } } //________________________________________________________________________ void AliAlignmentDataFilterITS::ConnectInputData(Option_t *) { // Connect ESD // Called once TTree* tree = dynamic_cast (GetInputData(0)); if(!tree) { printf("ERROR: Could not read chain from input slot 0\n"); } else { // Get the OCDB path and the list of OCDB objects used for reco TMap *cdbMap = (TMap*)(tree->GetTree()->GetUserInfo())->FindObject("cdbMap"); TList *cdbList = (TList*)(tree->GetTree()->GetUserInfo())->FindObject("cdbList"); //cdbList->Print(); // write the list to the user info of the output tree if(!fspTree) { printf("ERROR: fspTree does not exist\n"); } else { TMap *cdbMapCopy = new TMap(cdbMap->GetEntries()); cdbMapCopy->SetOwner(1); cdbMapCopy->SetName("cdbMap"); TIter iter1(cdbMap->GetTable()); TPair* pair = 0; while((pair = dynamic_cast (iter1.Next()))){ TObjString* keyStr = dynamic_cast (pair->Key()); TObjString* valStr = dynamic_cast (pair->Value()); cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName())); } TList *cdbListCopy = new TList(); cdbListCopy->SetOwner(1); cdbListCopy->SetName("cdbList"); TIter iter2(cdbList); TObjString* cdbEntry=0; while((cdbEntry =(TObjString*)(iter2.Next()))) { cdbListCopy->Add(new TObjString(*cdbEntry)); } cdbListCopy->Print(); fspTree->GetUserInfo()->Add(cdbMapCopy); fspTree->GetUserInfo()->Add(cdbListCopy); } // Disable all branches and enable only the needed ones tree->SetBranchStatus("fTriggerMask", 1); tree->SetBranchStatus("fSPDVertex*", 1); tree->SetBranchStatus("ESDfriend*", 1); tree->SetBranchAddress("ESDfriend.",&fESDfriend); tree->SetBranchStatus("fSPDMult*", 1); AliESDInputHandler *esdH = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); if(!esdH) { printf("ERROR: Could not get ESDInputHandler\n"); } else { fESD = esdH->GetEvent(); } } return; } //________________________________________________________________________ void AliAlignmentDataFilterITS::Init() { // Initialization return; } //________________________________________________________________________ void AliAlignmentDataFilterITS::CreateOutputObjects() { // Create the output container // // load the geometry if(!gGeoManager) { printf("AliAlignmentDataFilterITS::CreateOutputObjects(): loading geometry from %s\n",fGeometryFileName.Data()); AliGeomManager::LoadGeometry(fGeometryFileName.Data()); if(!gGeoManager) { printf("AliAlignmentDataFilterITS::CreateOutputObjects(): no geometry loaded \n"); return; } } // Several histograms are more conveniently managed in a TList fListOfHistos = new TList(); fListOfHistos->SetOwner(); fHistNevents = new TH1F("fHistNevents", "Number of processed events; N events; bin",5,-0.5,4.5); fHistNevents->Sumw2(); fHistNevents->SetMinimum(0); fListOfHistos->Add(fHistNevents); fHistNpoints = new TH1F("fHistNpoints", "Number of AliTrackPoints per track; N points; tracks",25,-0.5,24.5); fHistNpoints->Sumw2(); fHistNpoints->SetMinimum(0); fListOfHistos->Add(fHistNpoints); fHistPt = new TH1F("fHistPt", "p_{t} of tracks; p_{t} [GeV/c]; tracks",100,0,50); fHistPt->Sumw2(); fHistPt->SetMinimum(0); fListOfHistos->Add(fHistPt); Float_t zmax=14.; Int_t nbinsphi=20,nbinsz=4; fHistLayer0 = new TH2F("fHistLayer0","Points in layer inner SPD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); fListOfHistos->Add(fHistLayer0); zmax=14.; nbinsphi=40;nbinsz=4; fHistLayer1 = new TH2F("fHistLayer1","Points in layer outer SPD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); fListOfHistos->Add(fHistLayer1); zmax=22.; nbinsphi=14;nbinsz=6; fHistLayer2 = new TH2F("fHistLayer2","Points in layer inner SDD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); fListOfHistos->Add(fHistLayer2); zmax=29.5; nbinsphi=22;nbinsz=8; fHistLayer3 = new TH2F("fHistLayer3","Points in layer outer SDD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); fListOfHistos->Add(fHistLayer3); zmax=45.; nbinsphi=34;nbinsz=23; fHistLayer4 = new TH2F("fHistLayer4","Points in layer inner SSD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); fListOfHistos->Add(fHistLayer4); zmax=51.; nbinsphi=38;nbinsz=26; fHistLayer5 = new TH2F("fHistLayer5","Points in layer outer SSD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); fListOfHistos->Add(fHistLayer5); fntExtra = new TNtuple("fntExtra","extra clusters in ITS","ncls:layer:ladder:volid:phi:x:y:z:xloc:zloc:dxy:dz:d0mu:z0mu"); fListOfHistos->Add(fntExtra); fntCosmicMatching = new TNtuple("fntCosmicMatching","cosmic tracks matching in ITS","ncls1:ncls2:pt1:pt2:sigmad01:sigmad02:sigmaz01:sigmaz02:dxy:dz:phimu:thetamu:d0mu:z0mu"); fListOfHistos->Add(fntCosmicMatching); fspTree = new TTree("spTree","Tree with ITS track points"); const AliTrackPointArray *array = 0; Float_t curv,curverr,runNumber; const TObjString *itsaligndata = 0; fspTree->Branch("SP","AliTrackPointArray",&array); fspTree->Branch("curv",&curv); fspTree->Branch("curverr",&curverr); fspTree->Branch("run",&runNumber); fspTree->Branch("ITSAlignData",&itsaligndata); return; } //________________________________________________________________________ void AliAlignmentDataFilterITS::Exec(Option_t */*option*/) { // Execute analysis for current event: // write ITS AliTrackPoints for selected tracks to fspTree // check the geometry if(!gGeoManager) { printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n"); return; } // check if we have AliITSRecoParam if(!GetRecoParam()) { if(!fITSRecoParam) { printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n"); return; } } if(!fESD) { printf("AliAlignmentDataFilterITS::Exec(): no ESD \n"); return; } if(!fESDfriend) { printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n"); return; } // attach ESDfriend fESD->SetESDfriend(fESDfriend); // Post the data for slot 0 fHistNevents->Fill(0); // write field value to spTree UserInfo if(!((fspTree->GetUserInfo())->FindObject("BzkGauss"))) { Double_t bz=fESD->GetMagneticField(); TString bzString; bzString+=bz; TObjString *bzObjString = new TObjString(bzString); TList *bzList = new TList(); bzList->SetOwner(1); bzList->SetName("BzkGauss"); bzList->Add(bzObjString); fspTree->GetUserInfo()->Add(bzList); } // Process event as Cosmic or Collision //if(esd->GetEventType()== ???? ) { printf("AliAlignmentDataFilterITS::Exec(): MOVE ASAP TO esd->GetEventType() !\n"); if(GetRecoParam()->GetAlignFilterCosmics()) { FilterCosmic(fESD); } else { FilterCollision(fESD); } PostData(1,fListOfHistos); return; } //________________________________________________________________________ void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd) { // Extract ITS AliTrackPoints for Cosmics (check angular matching // of top and bottom track, merge the two tracks, if requested) // // Set branch addresses for space points tree AliTrackPointArray *arrayForTree=0; Float_t curv,curverr,runNumber; TObjString *itsaligndata=0; fspTree->SetBranchAddress("SP",&arrayForTree); fspTree->SetBranchAddress("curv",&curv); fspTree->SetBranchAddress("curverr",&curverr); fspTree->SetBranchAddress("run",&runNumber); fspTree->SetBranchAddress("ITSAlignData",&itsaligndata); runNumber = (Float_t)esd->GetRunNumber(); TTree* esdTree = dynamic_cast (GetInputData(0)); // Get the list of OCDB objects used for reco TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList"); TIter iter2(cdbList); TObjString* cdbEntry=0; TString cdbEntryString; while((cdbEntry =(TObjString*)(iter2.Next()))) { cdbEntryString = cdbEntry->GetString(); if(cdbEntryString.Contains("ITS/Align/Data")) itsaligndata = new TObjString(*cdbEntry); } TString triggeredClass = fESD->GetFiredTriggerClasses(); if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return; Int_t ntracks = esd->GetNumberOfTracks(); if(ntracks<2) return; if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return; Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos); Int_t *goodtracksArray = new Int_t[ntracks]; Float_t *phiArray = new Float_t[ntracks]; Float_t *thetaArray = new Float_t[ntracks]; Int_t *nclsArray = new Int_t[ntracks]; Int_t ngt=0; Int_t itrack=0; for (itrack=0; itrack < ntracks; itrack++) { AliESDtrack *track = esd->GetTrack(itrack); if (!track) continue; if(track->GetNcls(0)GetAlignFilterMinITSPoints()) continue; if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue; if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue; Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp()); Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl()); if(track->Pt()GetAlignFilterMinPt() || track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue; goodtracksArray[ngt] = itrack; phiArray[ngt] = phi; thetaArray[ngt] = theta; nclsArray[ngt] = track->GetNcls(0); ngt++; } if(ngt<2) { delete [] goodtracksArray; goodtracksArray=0; delete [] phiArray; phiArray=0; delete [] thetaArray; thetaArray=0; delete [] nclsArray; nclsArray=0; return; } // check matching of the two tracks from the muon Float_t min = 10000000.; Int_t maxCls = 0; Int_t good1 = -1, good2 = -1; for(Int_t itr1=0; itr1GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue; Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi()); if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue; if(nclsArray[itr1]+nclsArray[itr2] > maxCls) { maxCls = nclsArray[itr1]+nclsArray[itr2]; min = deltatheta+deltaphi; good1 = goodtracksArray[itr1]; good2 = goodtracksArray[itr2]; } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) { if(deltatheta+deltaphi < min) { min = deltatheta+deltaphi; good1 = goodtracksArray[itr1]; good2 = goodtracksArray[itr2]; } } } } delete [] goodtracksArray; goodtracksArray=0; delete [] phiArray; phiArray=0; delete [] thetaArray; thetaArray=0; delete [] nclsArray; nclsArray=0; if(good1<0) return; AliDebug(2,"ok track matching"); // track1 will be the inward track (top) // track2 the outward (bottom) AliESDtrack *track1=0; AliESDtrack *track2=0; AliESDtrack *track = esd->GetTrack(good1); if(track->Py()>0) { track1 = esd->GetTrack(good1); track2 = esd->GetTrack(good2); } else { track1 = esd->GetTrack(good2); track2 = esd->GetTrack(good1); } AliTrackPoint point; const AliTrackPointArray *array=0; Int_t ipt,volId,modId,layerId,lay,lad,det; Int_t jpt=0; Bool_t layerOK[6][2]; Int_t nclsTrk[2]={0,0}; for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE; for(itrack=0; itrack<2; itrack++) { if(itrack==0) { track = track1; } else { track = track2; } array = track->GetTrackPointArray(); if(!array) { AliWarning("No tracks points avaialble"); continue; } for(ipt=0; iptGetNPoints(); ipt++) { array->GetPoint(point,ipt); volId = point.GetVolumeID(); layerId = AliGeomManager::VolUIDToLayer(volId,modId); AliDebug(2,Form("%d %d\n",ipt,layerId-1)); if(point.IsExtra() && (GetRecoParam()->GetAlignFilterSkipExtra())) continue; if(layerId>6) continue; if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue; // check minAngleWrtITSModulePlanes Double_t p[3]; track->GetDirection(p); TVector3 pvec(p[0],p[1],p[2]); Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot); TVector3 normvec(rot[1],rot[4],rot[7]); Double_t angle = pvec.Angle(normvec); if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle; angle = 0.5*TMath::Pi()-angle; if(angleGetAlignFilterMinAngleWrtModulePlanes()) continue; layerOK[layerId-1][itrack]=kTRUE; jpt++; nclsTrk[itrack]++; } } AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1])); // read ITS cluster maps Int_t map1[6],map2[6]; for(Int_t ilay=0;ilay<6;ilay++) { map1[ilay]=0; map2[ilay]=0; if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1; if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1; } AliDebug(2,Form("ITS map 1: %d %d %d %d %d %d pt %f\n",map1[0],map1[1],map1[2],map1[3],map1[4],map1[5],track1->Pt())); AliDebug(2,Form("ITS map 2: %d %d %d %d %d %d pt %f\n",map2[0],map2[1],map2[2],map2[3],map2[4],map2[5],track2->Pt())); Int_t idx1[12],idx2[12]; track1->GetITSclusters(idx1); track2->GetITSclusters(idx2); AliDebug(2,Form("cls idx 1 %d %d %d %d %d %d %d %d %d %d %d %d\n",idx1[0],idx1[1],idx1[2],idx1[3],idx1[4],idx1[5],idx1[6],idx1[7],idx1[8],idx1[9],idx1[10],idx1[11])); AliDebug(2,Form("cls idx 2 %d %d %d %d %d %d %d %d %d %d %d %d\n",idx2[0],idx2[1],idx2[2],idx2[3],idx2[4],idx2[5],idx2[6],idx2[7],idx2[8],idx2[9],idx2[10],idx2[11])); if(jptGetAlignFilterMinITSPointsMerged()) return; AliDebug(2,Form(" Total points %d, accepted\n",jpt)); fHistNpoints->Fill(jpt); fHistPt->Fill(0.5*(track1->Pt()+track2->Pt())); Float_t d0z0mu[2]; track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu); //printf("d0mu %f z0mu %f\n",d0z0mu[0],d0z0mu[1]); Float_t dzOverlap[2]; Float_t curvArray[2],curverrArray[2]; Double_t globExtra[3],locExtra[3]; if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) arrayForTree = new AliTrackPointArray(jpt); jpt=0; for(itrack=0; itrack<2; itrack++) { if(itrack==0) { track = track1; } else { track = track2; } curvArray[itrack] = track->GetC(esd->GetMagneticField()); curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt(); if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) { jpt=0; arrayForTree = new AliTrackPointArray(nclsTrk[itrack]); } array = track->GetTrackPointArray(); for(ipt=0; iptGetNPoints(); ipt++) { array->GetPoint(point,ipt); volId = point.GetVolumeID(); layerId = AliGeomManager::VolUIDToLayer(volId,modId); if(layerId>6 || !layerOK[layerId-1][itrack]) continue; arrayForTree->AddPoint(jpt,&point); jpt++; switch(layerId) { case 1: fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 2: fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 3: fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 4: fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 5: fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 6: fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; } // Post the data for slot 0 if(jpt==1) PostData(1,fListOfHistos); // only if this is the first points if(!point.IsExtra() || !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue; nclsTrk[itrack]--; for(Int_t ll=1;llGetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap); AliTrackPoint pointT; Float_t radius,radiusT,phiv,phivT,thetav,thetavT; for(Int_t lll=0;lllGetPoint(pointT,lll); Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId); if(layerIdT!=layerId) continue; radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1])); radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1])); phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]); phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]); if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue; thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2])); thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2])); dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius)); if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue; dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT))); fntExtra->Fill((Float_t)nclsTrk[itrack],(Float_t)(layerId-1),lad,volId,TMath::ATan2(point.GetY(),point.GetX()),point.GetX(),point.GetY(),point.GetZ(),locExtra[0],locExtra[2],dzOverlap[0],dzOverlap[1],d0z0mu[0],d0z0mu[1]); } } if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) { curv = curvArray[itrack]; curverr = curverrArray[itrack]; fspTree->Fill(); } } if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) { curv = 0.5*(curvArray[0]+curvArray[1]); curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]); fspTree->Fill(); } PostData(0,fspTree); if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return; // fill ntuple with track-to-track matching Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty; Float_t d0[2],z0[2]; Float_t sigmad0[2],sigmaz0[2]; phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp()); thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl()); phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp()); thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl()); rotymu = TMath::ATan2(track1->Px(),track1->Pz()); rotyout = TMath::ATan2(track2->Px(),track2->Pz()); dphi = phimu - (phiout+TMath::Pi()); dtheta = thetamu - (TMath::Pi()-thetaout); if(rotymu>0) { droty = rotymu - (rotyout+TMath::Pi()); } else { droty = rotymu - (rotyout-TMath::Pi()); } Double_t alpha = TMath::ATan2(track1->Py(),track1->Px()); track1->Propagate(alpha,0.,esd->GetMagneticField()); track2->Propagate(alpha,0.,esd->GetMagneticField()); d0[0] = track1->GetY(); z0[0] = track1->GetZ(); d0[1] = track2->GetY(); z0[1] = track2->GetZ(); Float_t dxy = -(d0[0]-d0[1]); Float_t dz = z0[0]-z0[1]; sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2()); sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2()); sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2()); sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2()); /* Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3]; track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0); track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1); track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0); track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1); Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]); Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]); Float_t dxaty0 = x1aty0-x2aty0; */ fntCosmicMatching->Fill((Float_t)nclsTrk[0],(Float_t)nclsTrk[1],track1->Pt(),track2->Pt(),sigmad0[0],sigmad0[1],sigmaz0[0],sigmaz0[1],dxy,dz,phimu,thetamu,TMath::Abs(d0z0mu[0]),d0z0mu[1]); return; } //________________________________________________________________________ void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd) { // Extract ITS AliTrackPoints for Cosmics (check angular matching // of top and bottom track, merge the two tracks, if requested) // // Set branch addresses for space points tree AliTrackPointArray *arrayForTree=0; Float_t curv,curverr,runNumber; TObjString *itsaligndata=0; fspTree->SetBranchAddress("SP",&arrayForTree); fspTree->SetBranchAddress("curv",&curv); fspTree->SetBranchAddress("curverr",&curverr); fspTree->SetBranchAddress("run",&runNumber); fspTree->SetBranchAddress("ITSAlignData",&itsaligndata); runNumber = (Float_t)esd->GetRunNumber(); TTree* esdTree = dynamic_cast (GetInputData(0)); // Get the list of OCDB objects used for reco TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList"); TIter iter2(cdbList); TObjString* cdbEntry=0; TString cdbEntryString; while((cdbEntry =(TObjString*)(iter2.Next()))) { cdbEntryString = cdbEntry->GetString(); if(cdbEntryString.Contains("ITS/Align/Data")) itsaligndata = new TObjString(*cdbEntry); } Int_t ntracks = esd->GetNumberOfTracks(); if(ntracks==0) return; if(esd->GetPrimaryVertexTracks()->GetNContributors()<=0) return; Double_t vtxpos[3]; esd->GetPrimaryVertexTracks()->GetXYZ(vtxpos); Int_t ncls=0; Double_t pt=-10000.; Double_t d0z0[2],covd0z0[3]; const AliTrackPointArray *array = 0; for (Int_t itrack=0; itrack < ntracks; itrack++) { AliESDtrack * track = esd->GetTrack(itrack); if (!track) continue; if(track->GetNcls(0)GetAlignFilterMinITSPoints()) continue; if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue; if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue; if(track->Pt()GetAlignFilterMinPt() || track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue; pt = track->Pt(); ncls = track->GetNcls(0); Double_t maxd=10000.; track->PropagateToDCA(esd->GetPrimaryVertex(),esd->GetMagneticField(),maxd,d0z0,covd0z0); // read ITS cluster map Int_t map[6]; for(Int_t ilay=0;ilay<6;ilay++) { map[ilay]=0; if(track->HasPointOnITSLayer(ilay)) map[ilay]=1; } AliDebug(2,Form("ITS map : %d %d %d %d %d %d pt %f\n",map[0],map[1],map[2],map[3],map[4],map[5],track->Pt())); Int_t idx[12]; track->GetITSclusters(idx); AliDebug(2,Form("cls idx %d %d %d %d %d %d %d %d %d %d %d %d\n",idx[0],idx[1],idx[2],idx[3],idx[4],idx[5],idx[6],idx[7],idx[8],idx[9],idx[10],idx[11])); AliTrackPoint point; Int_t ipt,volId,modId,layerId,lay,lad,det; Int_t jpt=0; Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE; array = track->GetTrackPointArray(); if(!array) continue; for(ipt=0; iptGetNPoints(); ipt++) { array->GetPoint(point,ipt); volId = point.GetVolumeID(); layerId = AliGeomManager::VolUIDToLayer(volId,modId); if(layerId<1 || layerId>6) continue; if(point.IsExtra() && (GetRecoParam()->GetAlignFilterSkipExtra())) continue; layerOK[layerId-1]=kTRUE; jpt++; } if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue; fHistNpoints->Fill(jpt); fHistPt->Fill(pt); PostData(1,fListOfHistos); Float_t dzOverlap[2]; Double_t globExtra[3],locExtra[3]; arrayForTree = new AliTrackPointArray(jpt); jpt=0; array = track->GetTrackPointArray(); if(!array) continue; for(ipt=0; iptGetNPoints(); ipt++) { array->GetPoint(point,ipt); volId = point.GetVolumeID(); layerId = AliGeomManager::VolUIDToLayer(volId,modId); if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue; if(!point.IsExtra() || !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue; ncls--; for(Int_t ll=1;llGetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap); AliTrackPoint pointT; Float_t radius,radiusT,phiv,phivT,thetav,thetavT; for(Int_t lll=0;lllGetPoint(pointT,lll); Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId); if(layerIdT!=layerId) continue; radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1])); radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1])); phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]); phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]); if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue; thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2])); thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2])); dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius)); if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue; dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT))); fntExtra->Fill((Float_t)ncls,(Float_t)(layerId-1),lad,volId,TMath::ATan2(point.GetY(),point.GetX()),point.GetX(),point.GetY(),point.GetZ(),locExtra[0],locExtra[2],dzOverlap[0],dzOverlap[1],d0z0[0],d0z0[1]); } arrayForTree->AddPoint(jpt,&point); switch(layerId) { case 1: fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 2: fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 3: fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 4: fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 5: fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; case 6: fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); break; } jpt++; } curv = track->GetC(esd->GetMagneticField()); curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt(); fspTree->Fill(); } // end of tracks loop PostData(0,fspTree); return; } //________________________________________________________________________ void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/) { // Terminate analysis // AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n"); fspTree = dynamic_cast (GetOutputData(0)); if (!fspTree) { printf("ERROR: fspTree not available\n"); return; } fListOfHistos = dynamic_cast (GetOutputData(1)); if (!fListOfHistos) { printf("ERROR: fListOfHistos not available\n"); return; } fHistNevents = dynamic_cast(fListOfHistos->FindObject("fHistNevents")); fHistNpoints = dynamic_cast(fListOfHistos->FindObject("fHistNpoints")); fHistPt = dynamic_cast(fListOfHistos->FindObject("fHistPt")); fHistLayer0 = dynamic_cast(fListOfHistos->FindObject("fHistLayer0")); fHistLayer1 = dynamic_cast(fListOfHistos->FindObject("fHistLayer1")); fHistLayer2 = dynamic_cast(fListOfHistos->FindObject("fHistLayer2")); fHistLayer3 = dynamic_cast(fListOfHistos->FindObject("fHistLayer3")); fHistLayer4 = dynamic_cast(fListOfHistos->FindObject("fHistLayer4")); fHistLayer5 = dynamic_cast(fListOfHistos->FindObject("fHistLayer5")); fntExtra = dynamic_cast(fListOfHistos->FindObject("fntExtra")); fntCosmicMatching = dynamic_cast(fListOfHistos->FindObject("fntCosmicMatching")); return; }