/************************************************************************** * Copyright(c) 2007-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. * **************************************************************************/ //____________________________________________________________________ // // AliITSTrackleterSPDEff - find SPD chips efficiencies by using tracklets. // // This class has been developed from AliITSMultReconstructor (see // it for more details). It is the class for the Trackleter used to estimate // SPD plane efficiency. // The trackleter prediction is built using the vertex and 1 cluster. // // // Author : Giuseppe Eugenio Bruno, based on the skeleton of Reconstruct method provided by Tiziano Virgili // email: giuseppe.bruno@ba.infn.it // //____________________________________________________________________ /* $Id$ */ #include #include #include #include #include "AliITSMultReconstructor.h" #include "AliITSTrackleterSPDEff.h" #include "AliITSgeomTGeo.h" #include "AliLog.h" #include "AliITSPlaneEffSPD.h" #include "AliStack.h" //____________________________________________________________________ ClassImp(AliITSTrackleterSPDEff) //____________________________________________________________________ AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(): AliITSMultReconstructor(), fAssociationFlag1(0), fChipPredOnLay2(0), fChipPredOnLay1(0), fNTracklets1(0), fPhiWindowL1(0), fZetaWindowL1(0), fOnlyOneTrackletPerC1(0), fPlaneEffSPD(0), fMC(0), fUseOnlyPrimaryForPred(0), fUseOnlySecondaryForPred(0), fUseOnlySameParticle(0), fUseOnlyDifferentParticle(0), fUseOnlyStableParticle(0), fPredictionPrimary(0), fPredictionSecondary(0), fClusterPrimary(0), fClusterSecondary(0), fhClustersDPhiInterpAcc(0), fhClustersDThetaInterpAcc(0), fhClustersDZetaInterpAcc(0), fhClustersDPhiInterpAll(0), fhClustersDThetaInterpAll(0), fhClustersDZetaInterpAll(0), fhDPhiVsDThetaInterpAll(0), fhDPhiVsDThetaInterpAcc(0), fhDPhiVsDZetaInterpAll(0), fhDPhiVsDZetaInterpAcc(0), fhetaClustersLay2(0), fhphiClustersLay2(0) { // default constructor SetPhiWindowL1(); SetZetaWindowL1(); SetOnlyOneTrackletPerC1(); fAssociationFlag1 = new Bool_t[300000]; fChipPredOnLay2 = new UInt_t[300000]; fChipPredOnLay1 = new UInt_t[300000]; for(Int_t i=0; i<300000; i++) { fAssociationFlag1[i] = kFALSE; } if (GetHistOn()) BookHistos(); fPlaneEffSPD = new AliITSPlaneEffSPD(); } //______________________________________________________________________ AliITSTrackleterSPDEff::AliITSTrackleterSPDEff(const AliITSTrackleterSPDEff &mr) : AliITSMultReconstructor(mr), fAssociationFlag1(mr.fAssociationFlag1), fChipPredOnLay2(mr.fChipPredOnLay2), fChipPredOnLay1(mr.fChipPredOnLay1), fNTracklets1(mr.fNTracklets1), fPhiWindowL1(mr.fPhiWindowL1), fZetaWindowL1(mr.fZetaWindowL1), fOnlyOneTrackletPerC1(mr.fOnlyOneTrackletPerC1), fPlaneEffSPD(mr.fPlaneEffSPD), fMC(mr.fMC), fUseOnlyPrimaryForPred(mr.fUseOnlyPrimaryForPred), fUseOnlySecondaryForPred(mr.fUseOnlySecondaryForPred), fUseOnlySameParticle(mr.fUseOnlySameParticle), fUseOnlyDifferentParticle(mr.fUseOnlyDifferentParticle), fUseOnlyStableParticle(mr.fUseOnlyStableParticle), fPredictionPrimary(mr.fPredictionPrimary), fPredictionSecondary(mr.fPredictionSecondary), fClusterPrimary(mr.fClusterPrimary), fClusterSecondary(mr.fClusterSecondary), fhClustersDPhiInterpAcc(mr.fhClustersDPhiInterpAcc), fhClustersDThetaInterpAcc(mr.fhClustersDThetaInterpAcc), fhClustersDZetaInterpAcc(mr.fhClustersDZetaInterpAcc), fhClustersDPhiInterpAll(mr.fhClustersDPhiInterpAll), fhClustersDThetaInterpAll(mr.fhClustersDThetaInterpAll), fhClustersDZetaInterpAll(mr.fhClustersDZetaInterpAll), fhDPhiVsDThetaInterpAll(mr.fhDPhiVsDThetaInterpAll), fhDPhiVsDThetaInterpAcc(mr.fhDPhiVsDThetaInterpAcc), fhDPhiVsDZetaInterpAll(mr.fhDPhiVsDZetaInterpAll), fhDPhiVsDZetaInterpAcc(mr.fhDPhiVsDZetaInterpAcc), fhetaClustersLay2(mr.fhetaClustersLay2), fhphiClustersLay2(mr.fhphiClustersLay2) { // Copy constructor } //______________________________________________________________________ AliITSTrackleterSPDEff& AliITSTrackleterSPDEff::operator=(const AliITSTrackleterSPDEff& mr){ // Assignment operator this->~AliITSTrackleterSPDEff(); new(this) AliITSTrackleterSPDEff(mr); return *this; } //______________________________________________________________________ AliITSTrackleterSPDEff::~AliITSTrackleterSPDEff(){ // Destructor // delete histograms DeleteHistos(); delete [] fAssociationFlag1; delete [] fChipPredOnLay2; delete [] fChipPredOnLay1; delete [] fPredictionPrimary; delete [] fPredictionSecondary; delete [] fClusterPrimary; delete [] fClusterSecondary; // delete PlaneEff delete fPlaneEffSPD; } //____________________________________________________________________ void AliITSTrackleterSPDEff::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/,AliStack *pStack) { // // - calls LoadClusterArray that finds the position of the clusters // (in global coord) // - convert the cluster coordinates to theta, phi (seen from the // interaction vertex). Find the extrapolation/interpolation point. // - Find the chip corresponding to that // - Check if there is a cluster near that point // // reset counters fNClustersLay1 = 0; fNClustersLay2 = 0; fNTracklets = 0; fNSingleCluster = 0; // loading the clusters LoadClusterArrays(clusterTree); if(fMC && !pStack) {AliError("You asked for MC infos but AliStack not properly loaded"); return;} Bool_t found; Int_t nfTraPred1=0; Int_t ntTraPred1=0; Int_t nfTraPred2=0; Int_t ntTraPred2=0; Int_t nfClu1=0; Int_t ntClu1=0; Int_t nfClu2=0; Int_t ntClu2=0; // find the tracklets AliDebug(1,"Looking for tracklets... "); AliDebug(1,Form("Reconstruct: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2])); //########################################################### // Loop on layer 1 : finding theta, phi and z UInt_t key; for (Int_t iC1=0; iC1Fill(eta); fhphiClustersLay1->Fill(fClustersLay1[iC1][1]); } } // Loop on layer 2 : finding theta, phi and r for (Int_t iC2=0; iC2Fill(eta); fhphiClustersLay2->Fill(fClustersLay2[iC2][1]); } } //########################################################### // First part : Extrapolation to Layer 2 // Loop on layer 1 for (Int_t iC1=0; iC1TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; // find the difference in z (between linear projection from layer 1 // and the actual point: Dzeta= z1/r1*r2 -z2) Float_t r2 = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]); Float_t dZeta = TMath::Cos(fClustersLay1[iC1][0])*r2 - fClustersLay2[iC2][2]; if (fHistOn) { fhClustersDPhiAll->Fill(dPhi); fhClustersDThetaAll->Fill(dTheta); fhClustersDZetaAll->Fill(dZeta); fhDPhiVsDThetaAll->Fill(dTheta, dPhi); fhDPhiVsDZetaAll->Fill(dZeta, dPhi); } // make "elliptical" cut in Phi and Zeta! Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindow/fPhiWindow + dZeta*dZeta/fZetaWindow/fZetaWindow); if (d>1) continue; //look for the minimum distance: the minimum is in iC2WithBestDist if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) { distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi)); dPhimin = dPhi; dThetamin = dTheta; dZetamin = dZeta; iC2WithBestDist = iC2; } } } // end of loop over clusters in layer 2 if (distmin<100) { // This means that a cluster in layer 2 was found that matches with iC1 if (fHistOn) { fhClustersDPhiAcc->Fill(dPhimin); fhClustersDThetaAcc->Fill(dThetamin); fhClustersDZetaAcc->Fill(dZetamin); fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin); fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin); } if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association // store the tracklet // use the theta from the clusters in the first layer fTracklets[fNTracklets][0] = fClustersLay1[iC1][0]; // use the phi from the clusters in the first layer fTracklets[fNTracklets][1] = fClustersLay1[iC1][1]; // Store the difference between phi1 and phi2 fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1]; // find labels Int_t label1 = 0; Int_t label2 = 0; while (label2 < 3) { if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2]) break; label1++; if (label1 == 3) { label1 = 0; label2++; } } if (label2 < 3) { fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1]; } else { fTracklets[fNTracklets][3] = -2; } if (fHistOn) { Float_t eta=fTracklets[fNTracklets][0]; eta= TMath::Tan(eta/2.); eta=-TMath::Log(eta); fhetaTracklets->Fill(eta); fhphiTracklets->Fill(fTracklets[fNTracklets][1]); } // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation) found=FindChip(key,1,vtx,fClustersLay2[iC2WithBestDist][0],fClustersLay2[iC2WithBestDist][1],fClustersLay2[iC2WithBestDist][2]); if(!found){ AliWarning( Form("Reconstruct: cannot find chip on outer layer for cluster %d",iC2WithBestDist)); key=999999; } nfClu2+=(Int_t)found; // this for debugging purpose ntClu2++; // to check efficiency of the method FindChip if(key<1200) { // the Chip has been found if(fMC) { // this part only for MC // Int_t labc1=(Int_t)fClustersLay2[iC2WithBestDist][3]; // Int_t labc2=(Int_t)fClustersLay2[iC2WithBestDist][4]; // Int_t labc3=(Int_t)fClustersLay2[iC2WithBestDist][5]; if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it) if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it) } if (key==fChipPredOnLay2[iC1]) { // this control seems too loose: has to be checked ! // OK, success fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success } else { fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure // (might be in the tracking tollerance) } } fNTracklets++; } // if any cluster found --> increment statistics by 1 failure (provided you have chip prediction) else if (fChipPredOnLay2[iC1]<1200) fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay2[iC1]); } // end of loop over clusters in layer 1 fNTracklets1=fNTracklets; //################################################################### // Second part : Interpolation to Layer 1 // Loop on layer 2 for (Int_t iC2=0; iC2TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; // find the difference in z (between linear projection from layer 2 // and the actual point: Dzeta= z2/r2*r1 -z1) Float_t r1 = fClustersLay1[iC1][2]/TMath::Cos(fClustersLay1[iC1][0]); Float_t dZeta = TMath::Cos(fClustersLay2[iC2][0])*r1 - fClustersLay1[iC1][2]; if (fHistOn) { fhClustersDPhiInterpAll->Fill(dPhi); fhClustersDThetaInterpAll->Fill(dTheta); fhClustersDZetaInterpAll->Fill(dZeta); fhDPhiVsDThetaInterpAll->Fill(dTheta, dPhi); fhDPhiVsDZetaInterpAll->Fill(dZeta, dPhi); } // make "elliptical" cut in Phi and Zeta! Float_t d = TMath::Sqrt(dPhi*dPhi/fPhiWindowL1/fPhiWindowL1 + dZeta*dZeta/fZetaWindowL1/fZetaWindowL1); if (d>1) continue; //look for the minimum distance: the minimum is in iC1WithBestDist if (TMath::Sqrt(dZeta*dZeta+(r1*dPhi*r1*dPhi)) < distmin ) { distmin=TMath::Sqrt(dZeta*dZeta + (r1*dPhi*r1*dPhi)); dPhimin = dPhi; dThetamin = dTheta; dZetamin = dZeta; iC1WithBestDist = iC1; } } } // end of loop over clusters in layer 1 if (distmin<100) { // This means that a cluster in layer 1 was found that mathes with iC2 if (fHistOn) { fhClustersDPhiInterpAcc->Fill(dPhimin); fhClustersDThetaInterpAcc->Fill(dThetamin); fhClustersDZetaInterpAcc->Fill(dZetamin); fhDPhiVsDThetaInterpAcc->Fill(dThetamin, dPhimin); fhDPhiVsDZetaInterpAcc->Fill(dZetamin, dPhimin); } if (fOnlyOneTrackletPerC1) fAssociationFlag1[iC1WithBestDist] = kTRUE; // flag the association // flag the association // store the tracklet // use the theta from the clusters in the first layer fTracklets[fNTracklets][0] = fClustersLay2[iC2][0]; // use the phi from the clusters in the first layer fTracklets[fNTracklets][1] = fClustersLay2[iC2][1]; // Store the difference between phi1 and phi2 fTracklets[fNTracklets][2] = fClustersLay2[iC2][1] - fClustersLay1[iC1WithBestDist][1]; // find labels Int_t label1 = 0; Int_t label2 = 0; while (label2 < 3) { if ((Int_t) fClustersLay2[iC2][3+label1] != -2 && (Int_t) fClustersLay2[iC2][3+label1] == (Int_t) fClustersLay1[iC1WithBestDist][3+label2]) break; label1++; if (label1 == 3) { label1 = 0; label2++; } } if (label2 < 3) { fTracklets[fNTracklets][3] = fClustersLay2[iC2][3+label1]; } else { fTracklets[fNTracklets][3] = -2; } // Check that this cluster is still in the same chip (here you pass also Zvtx for better computation) found=FindChip(key,0,vtx,fClustersLay1[iC1WithBestDist][0],fClustersLay1[iC1WithBestDist][1],fClustersLay1[iC1WithBestDist][2]); if(!found){ AliWarning( Form("Reconstruct: cannot find chip on inner layer for cluster %d",iC1WithBestDist)); key=999999; } nfClu1+=(Int_t)found; // this for debugging purpose ntClu1++; // to check efficiency of the method FindChip if(key<1200) { if(fMC) { // this part only for MC // Int_t labc1=(Int_t)fClustersLay1[iC1WithBestDist][3]; // Int_t labc2=(Int_t)fClustersLay1[iC1WithBestDist][4]; // Int_t labc3=(Int_t)fClustersLay1[iC1WithBestDist][5]; if (fUseOnlyDifferentParticle && label2 < 3) continue; // same label (reject it) if (fUseOnlySameParticle && label2 == 3) continue; // different label (reject it) } if (key==fChipPredOnLay1[iC2]) { // this control seems too loose: has to be checked ! // OK, success fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // success } else { fPlaneEffSPD->UpDatePlaneEff(kTRUE,key); // this should not be a failure // (might be in the tracking tollerance) } } fNTracklets++; } // if no cluster found --> increment statistics by 1 failure (provided you have chip prediction) else if (fChipPredOnLay1[iC2]<1200) fPlaneEffSPD->UpDatePlaneEff(kFALSE,fChipPredOnLay1[iC2]); } // end of loop over clusters in layer 2 AliDebug(1,Form("%d tracklets found", fNTracklets)); AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 1 = %d / %d"),nfTraPred1,ntTraPred1)); AliDebug(1,Form(("Eff. of method FindChip for Track pred. on lay 2 = %d / %d"),nfTraPred2,ntTraPred2)); AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 1 = %d / %d"),nfClu1,ntClu1)); AliDebug(1,Form(("Eff. of method FindChip for Cluster on lay 2 = %d / %d"),nfClu2,ntClu2)); } //____________________________________________________________________ Bool_t AliITSTrackleterSPDEff::FindChip(UInt_t &key, Int_t layer, Float_t* vtx, Float_t thetaVtx, Float_t phiVtx, Float_t zVtx) { // // Input: a) layer number in the range [0,1] // b) vtx[3]: actual vertex // c) zVtx \ z of the cluster (-999 for tracklet) computed with respect to vtx // d) thetaVtx > theta and phi of the cluster/tracklet computed with respect to vtx // e) phiVtx / // Output: Unique key to locate a chip // return: kTRUE if succesfull if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return kFALSE;} Double_t r=GetRLayer(layer); //AliInfo(Form("Radius on layer %d is %f cm",layer,r)); // set phiVtx in the range [0,2pi] if(!SetAngleRange02Pi(phiVtx)) return kFALSE ; Double_t zAbs,phiAbs; // those are the polar coordinate, in the Absolute ALICE Reference // of the intersection of the tracklet with the pixel layer. if (TMath::Abs(zVtx)<100) zAbs=zVtx + vtx[2]; // this is fine only for the cluster, not for the track prediction else zAbs=r/TMath::Tan(thetaVtx) + vtx[2]; // this is the only way to do for the tracklet prediction AliDebug(1,Form("FindChip: vtx[0] = %f, vtx[1] = %f, vtx[2] = %f",vtx[0],vtx[1],vtx[2])); Double_t vtxy[2]={vtx[0],vtx[1]}; if (vtxy[0]*vtxy[1]+vtxy[1]*vtxy[1]>0) { // this method holds only for displaced vertices // this method gives you two interceptions if (!FindIntersectionPolar(vtxy,(Double_t)phiVtx,r,phiAbs)) return kFALSE; // set phiAbs in the range [0,2pi] if(!SetAngleRange02Pi(phiAbs)) return kFALSE; // since Vtx is very close to the ALICE origin, then phiVtx and phiAbs are very close; // therofore you can select the right intersection (among phiAbs1 and phiAbs2) by // taking the closest one to phiVtx AliDebug(1,Form("PhiVtx= %f, PhiAbs= %f",phiVtx,phiAbs)); } else phiAbs=phiVtx; Int_t idet=FindDetectorIndex(layer,phiAbs,zAbs); // this is the detector number // now you need to locate the chip within the idet detector, // starting from the local coordinates in such a detector Float_t locx; // local Cartesian coordinate (to be determined) corresponding to Float_t locz; // the Global Cilindrica coordinate (r,phiAbs,zAbs) . if(!FromGloCilToLocCart(layer,idet,r,phiAbs,zAbs, locx, locz)) return kFALSE; key=fPlaneEffSPD->GetKeyFromDetLocCoord(layer,idet,locx,locz); return kTRUE; } //______________________________________________________________________________ Double_t AliITSTrackleterSPDEff::GetRLayer(Int_t layer) { if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -999.;} Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 ! Double_t xyz[3], &x=xyz[0], &y=xyz[1]; AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz); Double_t r=TMath::Sqrt(x*x + y*y); AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz); r += TMath::Sqrt(x*x + y*y); AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz); r += TMath::Sqrt(x*x + y*y); AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz); r += TMath::Sqrt(x*x + y*y); r*=0.25; return r; } //______________________________________________________________________________ Bool_t AliITSTrackleterSPDEff::FromGloCilToLocCart(Int_t ilayer,Int_t idet, Double_t r, Double_t phi, Double_t z, Float_t &xloc, Float_t &zloc) { // this method transform Global Cilindrical coordinates into local (i.e. module) // cartesian coordinates // //Compute Cartesian Global Coordinate Double_t xyzGlob[3],xyzLoc[3]; xyzGlob[2]=z; xyzGlob[0]=r*TMath::Cos(phi); xyzGlob[1]=r*TMath::Sin(phi); xloc=0.; zloc=0.; if(idet<0) return kFALSE; Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6 Int_t lad = Int_t(idet/ndet) + 1; Int_t det = idet - (lad-1)*ndet + 1; AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc); xloc = (Float_t)xyzLoc[0]; zloc = (Float_t)xyzLoc[2]; return kTRUE; } //______________________________________________________________________________ Int_t AliITSTrackleterSPDEff::FindDetectorIndex(Int_t layer, Double_t phi, Double_t z) { //-------------------------------------------------------------------- //This function finds the detector crossed by the track //-------------------------------------------------------------------- if(layer<0 || layer >1) {AliWarning("Wrong layer: should be 0 or 1!"); return -1;} Int_t i=layer+1; // in AliITSgeomTGeo you count from 1 to 6 ! Int_t nladders=AliITSgeomTGeo::GetNLadders(i); Int_t ndetectors=AliITSgeomTGeo::GetNDetectors(i); Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z2=xyz[2]; AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz); Double_t phiOffset=TMath::ATan2(y,x); Double_t zOffset=z2; Double_t dphi; if (zOffset<0) // old geometry dphi = -(phi-phiOffset); else // new geometry dphi = phi-phiOffset; if (dphi < 0) dphi += 2*TMath::Pi(); else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi(); Int_t np=Int_t(dphi*nladders*0.5/TMath::Pi()+0.5); if (np>=nladders) np-=nladders; if (np<0) np+=nladders; Double_t dz=zOffset-z; Double_t nnz = dz*(ndetectors-1)*0.5/zOffset+0.5; Int_t nz = (nnz<0 ? -1 : (Int_t)nnz); if (nz>=ndetectors) {AliDebug(1,Form("too long: nz =%d",nz)); return -1;} if (nz<0) {AliDebug(1,Form("too short: nz =%d",nz)); return -1;} return np*ndetectors + nz; } //____________________________________________________________ Bool_t AliITSTrackleterSPDEff::FindIntersectionPolar(Double_t vtx[2],Double_t phiVtx, Double_t R,Double_t &phi) { // this method find the intersection in xy between a tracklet (straight line) and // a circonference (r=R), using polar coordinates. /* Input: - vtx[2]: actual vertex w.r.t. ALICE reference system - phiVtx: phi angle of the line (tracklet) computed w.r.t. vtx - R: radius of the circle Output: - phi : phi angle of the unique interception in the ALICE Global ref. system Correct method below: you have the equation of a circle (in polar coordinate) w.r.t. Actual vtx: r^2-2*r*r0*cos(phi-phi0) + r0^2 = R^2 , where (r0,phi0) is the centre of the circle In the same system, the equation of a semi-line is: phi=phiVtx; Hence you get one interception only: P=(r,phiVtx) Finally you want P in the ABSOLUTE ALICE system. */ Double_t rO=TMath::Sqrt(vtx[0]*vtx[0]+vtx[1]*vtx[1]); // polar coordinates of the ALICE origin Double_t phiO=TMath::ATan2(-vtx[1],-vtx[0]); // in the system with vtx[2] as origin Double_t bB=-2.*rO*TMath::Cos(phiVtx-phiO); Double_t cC=rO*rO-R*R; Double_t dDelta=bB*bB-4*cC; if(dDelta<0) return kFALSE; Double_t r1,r2; r1=(-bB-TMath::Sqrt(dDelta))/2; r2=(-bB+TMath::Sqrt(dDelta))/2; if(r1*r2>0) { printf("allora non hai capito nulla \n"); return kFALSE;} Double_t r=TMath::Max(r1,r2); // take the positive Double_t pvtx[2]; // Cartesian coordinates of the interception w.r.t. vtx Double_t pP[2]; // Cartesian coordinates of the interception w.r.t. ALICE origin pvtx[0]=r*TMath::Cos(phiVtx); pvtx[1]=r*TMath::Sin(phiVtx); pP[0]=vtx[0]+pvtx[0]; pP[1]=vtx[1]+pvtx[1]; phi=TMath::ATan2(pP[1],pP[0]); return kTRUE; } //___________________________________________________________ Bool_t AliITSTrackleterSPDEff::SetAngleRange02Pi(Double_t &angle) { while(angle >=2*TMath::Pi() || angle<0) { if(angle >= 2*TMath::Pi()) angle-=2*TMath::Pi(); if(angle < 0) angle+=2*TMath::Pi(); } return kTRUE; } //___________________________________________________________ Bool_t AliITSTrackleterSPDEff::PrimaryTrackChecker(Int_t ipart,AliStack* stack) { if(!fMC) {AliError("This method works only if SetMC() has been called"); return kFALSE;} if(!stack) {AliError("null pointer to MC stack"); return kFALSE;} if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return kFALSE;} // return stack->IsPhysicalPrimary(ipart); // looking at AliStack.cxx this does not seem to be complete (e.g. Pi0 Dalitz) if(!stack->IsPhysicalPrimary(ipart)) return kFALSE; // like below: as in the correction for Multiplicity (i.e. by hand in macro) TParticle* part = stack->Particle(ipart); TParticle* part0 = stack->Particle(0); // first primary TParticle* partl = stack->Particle(stack->GetNprimary()-1); //last primary if (part0->Vx()-partl->Vx()>0) AliDebug(1,Form("Difference in vtx position between 1th and last primaries %f %f %f", part0->Vx()-partl->Vx(),part0->Vy()-partl->Vy(), part0->Vz()-partl->Vz() )); if (!part || strcmp(part->GetName(),"XXX")==0) {AliWarning("String , not particle ??") ;return kFALSE; } TParticlePDG* pdgPart = part->GetPDG(); if (TMath::Abs(pdgPart->Charge()) < 3) {AliWarning("This seems a quark"); return kFALSE;} Double_t distx = part->Vx() - part0->Vx(); Double_t disty = part->Vy() - part0->Vy(); Double_t distz = part->Vz() - part0->Vz(); Double_t distR=TMath::Sqrt(distx*distx + disty*disty + distz*distz); if (distR > 0.05) {AliDebug(1,Form("True vertex should be %f %f, this particle from %f %f ", part0->Vx(),part0->Vy(),part->Vx(),part->Vy())); return kFALSE; }// primary if within 500 microns from true Vertex if(fUseOnlyStableParticle && DecayingTrackChecker(ipart,stack)<2) return kFALSE; return kTRUE; } //_____________________________________________________________________________________________ Int_t AliITSTrackleterSPDEff::DecayingTrackChecker(Int_t ipart,AliStack* stack) { if(!fMC) {AliError("This method works only if SetMC() has been called"); return 0;} if(!stack) {AliError("null pointer to MC stack"); return 0;} if(ipart >= stack->GetNtrack()) {AliError("this track label is not in MC stack"); return 0;} TParticle* part = stack->Particle(ipart); //TParticle* part0 = stack->Particle(0); // first primary Int_t nret=0; TParticle* dau = 0; Int_t nDau = 0; Int_t firstDau = part->GetFirstDaughter(); if (firstDau > 0) { Int_t lastDau = part->GetLastDaughter(); nDau = lastDau - firstDau + 1; //printf("number of daugthers %d \n",nDau); if (nDau > 0) { //for(Int_t j=firstDau; j<=lastDau; j++) for(Int_t j=firstDau; j<=firstDau; j++) { // only first one dau = stack->Particle(j); Double_t distx = dau->Vx()-part->Vx(); Double_t disty = dau->Vy()-part->Vy(); Double_t distz = dau->Vz()-part->Vz(); Double_t distR = TMath::Sqrt(distx*distx+disty*disty+distz*distz); if (distR > GetRLayer(0)+0.5) nret=1; // decay after first pixel layer if (distR > GetRLayer(1)+0.5) nret=2; // decay after second pixel layer } } } else nret = 3; // stable particle return nret; } //_________________________________________________________________ void AliITSTrackleterSPDEff::InitPredictionMC() { if(!fMC) {AliError("This method works only if SetMC() has been called"); return;} fPredictionPrimary = new Int_t[1200]; fPredictionSecondary = new Int_t[1200]; fClusterPrimary = new Int_t[1200]; fClusterSecondary = new Int_t[1200]; for(Int_t i=0; i<1200; i++) { fPredictionPrimary[i]=0; fPredictionSecondary[i]=0; fPredictionSecondary[i]=0; fClusterSecondary[i]=0; } return; } //______________________________________________________________________ Int_t AliITSTrackleterSPDEff::GetPredictionPrimary(const UInt_t key) const { // // This method return the Data menmber fPredictionPrimary [1200]. // You can call it only for MC events. // fPredictionPrimary[key] contains the number of tracklet predictions on the // given chip key built using a cluster on the other layer produced (at least) // from a primary particle. // Key refers to the chip crossed by the prediction // // if (!fMC) {CallWarningMC(); return 0;} if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;} return fPredictionPrimary[(Int_t)key]; } //______________________________________________________________________ Int_t AliITSTrackleterSPDEff::GetPredictionSecondary(const UInt_t key) const { // // This method return the Data menmber fPredictionSecondary [1200]. // You can call it only for MC events. // fPredictionSecondary[key] contains the number of tracklet predictions on the // given chip key built using a cluster on the other layer produced (only) // from a secondary particle // Key refers to the chip crossed by the prediction // // if (!fMC) {CallWarningMC(); return 0;} if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;} return fPredictionSecondary[(Int_t)key]; } //______________________________________________________________________ Int_t AliITSTrackleterSPDEff::GetClusterPrimary(const UInt_t key) const { // // This method return the Data menmber fClusterPrimary [1200]. // You can call it only for MC events. // fClusterPrimary[key] contains the number of tracklet predictions // built using a cluster on that layer produced (only) // from a primary particle // Key refers to the chip used to build the prediction // // if (!fMC) {CallWarningMC(); return 0;} if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;} return fClusterPrimary[(Int_t)key]; } //______________________________________________________________________ Int_t AliITSTrackleterSPDEff::GetClusterSecondary(const UInt_t key) const { // // This method return the Data menmber fClusterSecondary [1200]. // You can call it only for MC events. // fClusterSecondary[key] contains the number of tracklet predictions // built using a cluster on that layer produced (only) // from a secondary particle // Key refers to the chip used to build the prediction // if (!fMC) {CallWarningMC(); return 0;} if (key>=1200) {AliWarning("You asked for a non existing chip"); return -999;} return fClusterSecondary[(Int_t)key]; } //______________________________________________________________________ void AliITSTrackleterSPDEff::PrintAscii(ostream *os)const{ // Print out some class data values in Ascii Form to output stream // Inputs: // ostream *os Output stream where Ascii data is to be writen // Outputs: // none. // Return: // none. *os << fPhiWindowL1 <<" "<< fZetaWindowL1 << " " << fPhiWindow <<" "<< fZetaWindow ; *os << " " << fMC; if(!fMC) {AliInfo("Writing only cuts, no MC info"); return;} *os << " " << fUseOnlyPrimaryForPred << " " << fUseOnlySecondaryForPred << " " << fUseOnlySameParticle << " " << fUseOnlyDifferentParticle << " " << fUseOnlyStableParticle ; for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionPrimary(i) ; for(Int_t i=0;i<1200;i++) *os <<" "<< GetPredictionSecondary(i) ; for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterPrimary(i) ; for(Int_t i=0;i<1200;i++) *os <<" "<< GetClusterSecondary(i) ; return; } //______________________________________________________________________ void AliITSTrackleterSPDEff::ReadAscii(istream *is){ // Read in some class data values in Ascii Form to output stream // Inputs: // istream *is Input stream where Ascii data is to be read in from // Outputs: // none. // Return: // none. *is >> fPhiWindowL1 >> fZetaWindowL1 >> fPhiWindow >> fZetaWindow; *is >> fMC; if(!fMC) {AliInfo("Reading only cuts, no MC info available");return;} *is >> fUseOnlyPrimaryForPred >> fUseOnlySecondaryForPred >> fUseOnlySameParticle >> fUseOnlyDifferentParticle >> fUseOnlyStableParticle; for(Int_t i=0;i<1200;i++) *is >> fPredictionPrimary[i] ; for(Int_t i=0;i<1200;i++) *is >> fPredictionSecondary[i] ; for(Int_t i=0;i<1200;i++) *is >> fClusterPrimary[i] ; for(Int_t i=0;i<1200;i++) *is >> fClusterSecondary[i] ; return; } //______________________________________________________________________ ostream &operator<<(ostream &os,const AliITSTrackleterSPDEff &s){ // Standard output streaming function // Inputs: // ostream &os output steam // AliITSTrackleterSPDEff &s class to be streamed. // Output: // none. // Return: // ostream &os The stream pointer s.PrintAscii(&os); return os; } //______________________________________________________________________ istream &operator>>(istream &is,AliITSTrackleterSPDEff &s){ // Standard inputput streaming function // Inputs: // istream &is input steam // AliITSTrackleterSPDEff &s class to be streamed. // Output: // none. // Return: // ostream &os The stream pointer //printf("prova %d \n", (Int_t)s.GetMC()); s.ReadAscii(&is); return is; } //______________________________________________________________________ void AliITSTrackleterSPDEff::SavePredictionMC(TString filename) const { // // This Method write into an asci file (do not know why binary does not work) // the used cuts and the statistics of the MC related quantities // The method SetMC() has to be called before // Input TString filename: name of file for output (it deletes already existing // file) // Output: none // // if(!fMC) {CallWarningMC(); return;} ofstream out(filename.Data(),ios::out | ios::binary); out << *this; out.close(); return; } //____________________________________________________________________ void AliITSTrackleterSPDEff::ReadPredictionMC(TString filename) { // // This Method read from an asci file (do not know why binary does not work) // the cuts to be used and the statistics of the MC related quantities // Input TString filename: name of input file for output // The method SetMC() has to be called before // Output: none // // if(!fMC) {CallWarningMC(); return;} if( gSystem->AccessPathName( filename.Data() ) ) { AliError( Form( "file (%s) not found", filename.Data() ) ); return; } ifstream in(filename.Data(),ios::in | ios::binary); in >> *this; in.close(); return; } //____________________________________________________________________ Bool_t AliITSTrackleterSPDEff::SaveHists() { // This (private) method save the histograms on the output file // (only if fHistOn is TRUE). // Also the histograms from the base class are saved through the // AliITSMultReconstructor::SaveHists() call if (!GetHistOn()) return kFALSE; AliITSMultReconstructor::SaveHists(); // this save the histograms of the base class fhClustersDPhiInterpAll->Write(); fhClustersDThetaInterpAll->Write(); fhClustersDZetaInterpAll->Write(); fhDPhiVsDThetaInterpAll->Write(); fhDPhiVsDZetaInterpAll->Write(); fhClustersDPhiInterpAcc->Write(); fhClustersDThetaInterpAcc->Write(); fhClustersDZetaInterpAcc->Write(); fhDPhiVsDThetaInterpAcc->Write(); fhDPhiVsDZetaInterpAcc->Write(); fhetaClustersLay2->Write(); fhphiClustersLay2->Write(); return kTRUE; } //__________________________________________________________ Bool_t AliITSTrackleterSPDEff::WriteHistosToFile(TString filename, Option_t* option) { // // Saves the histograms into a tree and saves the trees into a file // Also the histograms from the base class are saved // if (!GetHistOn()) return kFALSE; if (filename.Data()=="") { AliWarning("WriteHistosToFile: null output filename!"); return kFALSE; } TFile *hFile=new TFile(filename.Data(),option, "The File containing the histos for SPD efficiency studies with tracklets"); if(!SaveHists()) return kFALSE; hFile->Write(); hFile->Close(); return kTRUE; } //____________________________________________________________ void AliITSTrackleterSPDEff::BookHistos() { // // This method books addtitional histograms // w.r.t. those of the base class. // In particular, the differences of cluster coordinate between the two SPD // layers are computed in the interpolation phase // if (! GetHistOn()) { AliInfo("Call SetHistOn(kTRUE) first"); return;} fhClustersDPhiInterpAcc = new TH1F("dphiaccInterp", "dphi Interpolation phase", 100,0.,0.1); fhClustersDPhiInterpAcc->SetDirectory(0); fhClustersDThetaInterpAcc = new TH1F("dthetaaccInterp","dtheta Interpolation phase",100,-0.1,0.1); fhClustersDThetaInterpAcc->SetDirectory(0); fhClustersDZetaInterpAcc = new TH1F("dzetaaccInterp","dzeta Interpolation phase",100,-1.,1.); fhClustersDZetaInterpAcc->SetDirectory(0); fhDPhiVsDZetaInterpAcc = new TH2F("dphiVsDzetaaccInterp","dphiVsDzeta Interpolation phase",100,-1.,1.,100,0.,0.1); fhDPhiVsDZetaInterpAcc->SetDirectory(0); fhDPhiVsDThetaInterpAcc = new TH2F("dphiVsDthetaAccInterp","dphiVsDtheta Interpolation phase",100,-0.1,0.1,100,0.,0.1); fhDPhiVsDThetaInterpAcc->SetDirectory(0); fhClustersDPhiInterpAll = new TH1F("dphiallInterp", "dphi Interpolation phase", 100,0.0,0.5); fhClustersDPhiInterpAll->SetDirectory(0); fhClustersDThetaInterpAll = new TH1F("dthetaallInterp","dtheta Interpolation phase",100,-0.5,0.5); fhClustersDThetaInterpAll->SetDirectory(0); fhClustersDZetaInterpAll = new TH1F("dzetaallInterp","dzeta Interpolation phase",100,-5.,5.); fhClustersDZetaInterpAll->SetDirectory(0); fhDPhiVsDZetaInterpAll = new TH2F("dphiVsDzetaallInterp","dphiVsDzeta Interpolation phase",100,-5.,5.,100,0.,0.5); fhDPhiVsDZetaInterpAll->SetDirectory(0); fhDPhiVsDThetaInterpAll = new TH2F("dphiVsDthetaAllInterp","dphiVsDtheta Interpolation phase",100,-0.5,0.5,100,0.,0.5); fhDPhiVsDThetaInterpAll->SetDirectory(0); fhetaClustersLay2 = new TH1F("etaClustersLay2", "etaCl2", 100,-2.,2.); fhetaClustersLay2->SetDirectory(0); fhphiClustersLay2 = new TH1F("phiClustersLay2", "phiCl2", 100, 0., 2*TMath::Pi()); fhphiClustersLay2->SetDirectory(0); return; } //____________________________________________________________ void AliITSTrackleterSPDEff::DeleteHistos() { // // Private method to delete Histograms from memory // it is called. e.g., by the destructor. // if(fhClustersDPhiInterpAcc) {delete fhClustersDPhiInterpAcc; fhClustersDPhiInterpAcc=0;} if(fhClustersDThetaInterpAcc) {delete fhClustersDThetaInterpAcc; fhClustersDThetaInterpAcc=0;} if(fhClustersDZetaInterpAcc) {delete fhClustersDZetaInterpAcc; fhClustersDZetaInterpAcc=0;} if(fhClustersDPhiInterpAll) {delete fhClustersDPhiInterpAll; fhClustersDPhiInterpAll=0;} if(fhClustersDThetaInterpAll) {delete fhClustersDThetaInterpAll; fhClustersDThetaInterpAll=0;} if(fhClustersDZetaInterpAll) {delete fhClustersDZetaInterpAll; fhClustersDZetaInterpAll=0;} if(fhDPhiVsDThetaInterpAll) {delete fhDPhiVsDThetaInterpAll; fhDPhiVsDThetaInterpAll=0;} if(fhDPhiVsDThetaInterpAcc) {delete fhDPhiVsDThetaInterpAcc; fhDPhiVsDThetaInterpAcc=0;} if(fhDPhiVsDZetaInterpAll) {delete fhDPhiVsDZetaInterpAll; fhDPhiVsDZetaInterpAll=0;} if(fhDPhiVsDZetaInterpAcc) {delete fhDPhiVsDZetaInterpAcc; fhDPhiVsDZetaInterpAcc=0;} if(fhetaClustersLay2) {delete fhetaClustersLay2; fhetaClustersLay2=0;} if(fhphiClustersLay2) {delete fhphiClustersLay2; fhphiClustersLay2=0;} } //_______________________________________________________________