+
+ return;
+}
+
+//_____________________________________________________________________________
+void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
+ Int_t minITSpts,Float_t maxMatchingAngle,
+ Bool_t cuts,
+ Float_t minAngleWrtITSModulePlanes,
+ Float_t minMom,Float_t maxMom,
+ Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
+ Float_t minSinTheta,Float_t maxSinTheta)
+{
+ // Analyzes and filters ESD tracks
+ // Merges inward and outward tracks in one single track
+ // Stores the selected track space points
+ // into the output file
+
+ if (!fESDChain) return;
+
+ AliESDEvent *esd = new AliESDEvent();
+ esd->ReadFromTree(fESDChain);
+ AliESDfriend *esdf = 0;
+ fESDChain->SetBranchStatus("ESDfriend*",1);
+ fESDChain->SetBranchAddress("ESDfriend.",&esdf);
+
+ // Open the output file
+ if (fPointsFilename.IsNull()) {
+ AliWarning("Incorrect output filename!");
+ return;
+ }
+
+ TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
+ if (!pointsFile || !pointsFile->IsOpen()) {
+ AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
+ return;
+ }
+
+ TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
+ const AliTrackPointArray *array = 0;
+ AliTrackPointArray *array2 = 0;
+ pointsTree->Branch("SP","AliTrackPointArray", &array2);
+
+ Int_t ievent = 0;
+ while (fESDChain->GetEntry(ievent++)) {
+ if (!esd) break;
+
+ esd->SetESDfriend(esdf); //Attach the friend to the ESD
+
+ Int_t ntracks = esd->GetNumberOfTracks();
+ if(ntracks<2) continue;
+ Int_t *goodtracksArray = new Int_t[ntracks];
+ Float_t *phiArray = new Float_t[ntracks];
+ Float_t *thetaArray = new Float_t[ntracks];
+ Int_t ngt=0;
+ for (Int_t itrack=0; itrack < ntracks; itrack++) {
+ AliESDtrack * track = esd->GetTrack(itrack);
+ if (!track) continue;
+
+ if(track->GetNcls(0) < minITSpts) continue;
+ Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
+ Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
+ if(cuts) {
+ if(track->GetP()<minMom || track->GetP()>maxMom) continue;
+ Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
+ if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
+ Float_t sintheta = TMath::Sin(theta);
+ if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
+ }
+ goodtracksArray[ngt]=itrack;
+ phiArray[ngt]=phi;
+ thetaArray[ngt]=theta;
+ ngt++;
+ }
+
+ if(ngt<2) {
+ delete [] goodtracksArray; goodtracksArray=0;
+ delete [] phiArray; phiArray=0;
+ delete [] thetaArray; thetaArray=0;
+ continue;
+ }
+
+ // check matching of the two tracks from the muon
+ Float_t min = 10000000.;
+ Int_t good1 = -1, good2 = -1;
+ for(Int_t itr1=0; itr1<ngt-1; itr1++) {
+ for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
+ Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
+ if(deltatheta>maxMatchingAngle) continue;
+ Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
+ if(deltaphi>maxMatchingAngle) continue;
+ //printf("%f %f %f %f\n",deltaphi,deltatheta,thetaArray[itr1],thetaArray[itr2]);
+ 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;
+
+ if(good1<0) continue;
+
+ AliESDtrack * track1 = esd->GetTrack(good1);
+ AliESDtrack * track2 = esd->GetTrack(good2);
+
+ AliTrackPoint point;
+ Int_t ipt,volId,modId,layerId;
+ Int_t jpt=0;
+ Bool_t layerOK[6][2];
+ for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kTRUE;
+ array = track1->GetTrackPointArray();
+ for(ipt=0; ipt<array->GetNPoints(); ipt++) {
+ array->GetPoint(point,ipt);
+ if(onlyITS) {
+ volId = point.GetVolumeID();
+ layerId = AliGeomManager::VolUIDToLayer(volId,modId);
+ if(layerId>6) continue;
+ // check minAngleWrtITSModulePlanes
+ if(cuts) {
+ Double_t p[3]; track1->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(angle<minAngleWrtITSModulePlanes) {
+ layerOK[layerId-1][0]=kFALSE;
+ continue;
+ }
+ }
+ }
+ jpt++;
+ }
+ array = track2->GetTrackPointArray();
+ for(ipt=0; ipt<array->GetNPoints(); ipt++) {
+ array->GetPoint(point,ipt);
+ if(onlyITS) {
+ volId = point.GetVolumeID();
+ layerId = AliGeomManager::VolUIDToLayer(volId,modId);
+ if(layerId>6) continue;
+ // check minAngleWrtITSModulePlanes
+ if(cuts) {
+ Double_t p[3]; track2->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(angle<minAngleWrtITSModulePlanes) {
+ layerOK[layerId-1][0]=kFALSE;
+ continue;
+ }
+ }
+ }
+ jpt++;
+ }
+
+ if(jpt < 2*minITSpts) continue;
+ array2 = new AliTrackPointArray(jpt);
+ jpt=0;
+ array = track1->GetTrackPointArray();
+ for(ipt=0; ipt<array->GetNPoints(); ipt++) {
+ array->GetPoint(point,ipt);
+ if(onlyITS) {
+ volId = point.GetVolumeID();
+ layerId = AliGeomManager::VolUIDToLayer(volId,modId);
+ if(layerId>6 || !layerOK[layerId-1][0]) continue;
+ }
+ array2->AddPoint(jpt,&point);
+ jpt++;
+ }
+ array = track2->GetTrackPointArray();
+ for(ipt=0; ipt<array->GetNPoints(); ipt++) {
+ array->GetPoint(point,ipt);
+ if(onlyITS) {
+ volId = point.GetVolumeID();
+ layerId = AliGeomManager::VolUIDToLayer(volId,modId);
+ if(layerId>6 || !layerOK[layerId-1][1]) continue;
+ }
+ array2->AddPoint(jpt,&point);
+ jpt++;
+ }
+
+ pointsTree->Fill();
+ }
+
+ if (!pointsTree->Write()) {
+ AliWarning("Can't write the tree with track point arrays!");
+ return;
+ }
+
+ pointsFile->Close();
+ return;