PWG1/AliAlignmentDataFilterITS.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTask to extract from ESD tracks the AliTrackPointArrays
  19 // with ITS points for selected tracks. This are the input data for alignment
  20 //
  21 // Author: A.Dainese, andrea.dainese@pd.infn.it
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TTree.h>
  25 #include <TChain.h>
  26 #include <TNtuple.h>
  27 #include <TList.h>
  28 #include <TH1F.h>
  29 #include <TH2F.h>
  30 #include <TVector3.h>
  31 #include <TGeoManager.h>
  32
  33 #include "AliLog.h"
  34 #include "AliGeomManager.h"
  35 #include "AliITSgeomTGeo.h"
  36 #include "AliTrackPointArray.h"
  37 #include "AliESDInputHandler.h"
  38 #include "AliESDVertex.h"
  39 #include "AliESDtrack.h"
  40 #include "AliESDEvent.h"
  41 #include "AliESDfriend.h"
  42 #include "AliAnalysisTask.h"
  43 #include "AliAnalysisManager.h"
  44 #include "AliAlignmentDataFilterITS.h"
  45
  46
  47 ClassImp(AliAlignmentDataFilterITS)
  48
  49
  50 //________________________________________________________________________
  51 AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name):
  52 AliAnalysisTask(name,"task"),
  53 fOnlySPDFO(kFALSE),
  54 fGeometryFileName("geometry.root"),
  55 fITSRecoParam(0),
  56 fESD(0),
  57 fESDfriend(0),
  58 fListOfHistos(0),
  59 fspTree(0),
  60 fHistNpoints(0),
  61 fHistPt(0),
  62 fHistLayer0(0),
  63 fHistLayer1(0),
  64 fHistLayer2(0),
  65 fHistLayer3(0),
  66 fHistLayer4(0),
  67 fHistLayer5(0),
  68 fntExtra(0),
  69 fntCosmicMatching(0)
  70 {
  71   // Constructor
  72
  73   // Define input and output slots here
  74   // Input slot #0 works with a TChain
  75   DefineInput(0, TChain::Class());
  76   // Output slot #0 writes into a TTree
  77   DefineOutput(0,TTree::Class());  //My private output
  78   // Output slot #1 writes into a TList
  79   DefineOutput(1,TList::Class());  //My private output
  80 }
  81
  82 //________________________________________________________________________
  83 AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS()
  84 {
  85   // Destructor
  86   if (fListOfHistos) {
  87     delete fListOfHistos;
  88     fListOfHistos = 0;
  89   }
  90   if (fspTree) {
  91     delete fspTree;
  92     fspTree = 0;
  93   }
  94   if (fHistNpoints) {
  95     delete fHistNpoints;
  96     fHistNpoints = 0;
  97   }
  98   if (fHistPt) {
  99     delete fHistPt;
 100     fHistPt = 0;
 101   }
 102   if (fHistLayer0) {
 103     delete fHistLayer0;
 104     fHistLayer0 = 0;
 105   }
 106   if (fHistLayer1) {
 107     delete fHistLayer1;
 108     fHistLayer1 = 0;
 109   }
 110   if (fHistLayer2) {
 111     delete fHistLayer2;
 112     fHistLayer2 = 0;
 113   }
 114   if (fHistLayer3) {
 115     delete fHistLayer3;
 116     fHistLayer3 = 0;
 117   }
 118   if (fHistLayer4) {
 119     delete fHistLayer4;
 120     fHistLayer4 = 0;
 121   }
 122   if (fHistLayer5) {
 123     delete fHistLayer5;
 124     fHistLayer5 = 0;
 125   }
 126   if (fntExtra) {
 127     delete fntExtra;
 128     fntExtra = 0;
 129   }
 130   if (fntCosmicMatching) {
 131     delete fntCosmicMatching;
 132     fntCosmicMatching = 0;
 133   }
 134 }
 135 //________________________________________________________________________
 136 void AliAlignmentDataFilterITS::ConnectInputData(Option_t *)
 137 {
 138   // Connect ESD
 139   // Called once
 140
 141   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 142   if(!tree) {
 143     printf("ERROR: Could not read chain from input slot 0\n");
 144   } else {
 145     // Disable all branches and enable only the needed ones
 146
 147     tree->SetBranchStatus("fTriggerMask", 1);
 148     tree->SetBranchStatus("fSPDVertex*", 1);
 149
 150     tree->SetBranchStatus("ESDfriend*", 1);
 151     tree->SetBranchAddress("ESDfriend.",&fESDfriend);
 152
 153     tree->SetBranchStatus("fSPDMult*", 1);
 154
 155     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 156
 157     if(!esdH) {
 158       printf("ERROR: Could not get ESDInputHandler\n");
 159     } else {
 160       fESD = esdH->GetEvent();
 161     }
 162   }
 163
 164   return;
 165 }
 166
 167 //________________________________________________________________________
 168 void AliAlignmentDataFilterITS::Init()
 169 {
 170   // Initialization
 171
 172   return;
 173 }
 174
 175 //________________________________________________________________________
 176 void AliAlignmentDataFilterITS::CreateOutputObjects()
 177 {
 178   // Create the output container
 179   //
 180
 181   // Several histograms are more conveniently managed in a TList
 182   fListOfHistos = new TList();
 183   fListOfHistos->SetOwner();
 184
 185   fHistNpoints = new TH1F("fHistNpoints", "Number of AliTrackPoints per track; N points; tracks",25,-0.5,24.5);
 186   fHistNpoints->Sumw2();
 187   fHistNpoints->SetMinimum(0);
 188   fListOfHistos->Add(fHistNpoints);
 189
 190   fHistPt = new TH1F("fHistPt", "p_{t} of tracks; p_{t} [GeV/c]; tracks",100,0,50);
 191   fHistPt->Sumw2();
 192   fHistPt->SetMinimum(0);
 193   fListOfHistos->Add(fHistPt);
 194
 195
 196   Float_t zmax=14.;
 197   Int_t nbinsphi=20,nbinsz=4;
 198   fHistLayer0 = new TH2F("fHistLayer0","Points in layer inner SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 199   fListOfHistos->Add(fHistLayer0);
 200   zmax=14.;
 201   nbinsphi=40;nbinsz=4;
 202   fHistLayer1 = new TH2F("fHistLayer1","Points in layer outer SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 203   fListOfHistos->Add(fHistLayer1);
 204   zmax=22.;
 205   nbinsphi=14;nbinsz=6;
 206   fHistLayer2 = new TH2F("fHistLayer2","Points in layer inner SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 207   fListOfHistos->Add(fHistLayer2);
 208   zmax=29.5;
 209   nbinsphi=22;nbinsz=8;
 210   fHistLayer3 = new TH2F("fHistLayer3","Points in layer outer SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 211   fListOfHistos->Add(fHistLayer3);
 212   zmax=45.;
 213   nbinsphi=34;nbinsz=23;
 214   fHistLayer4 = new TH2F("fHistLayer4","Points in layer inner SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 215   fListOfHistos->Add(fHistLayer4);
 216   zmax=51.;
 217   nbinsphi=38;nbinsz=26;
 218   fHistLayer5 = new TH2F("fHistLayer5","Points in layer outer SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 219   fListOfHistos->Add(fHistLayer5);
 220
 221
 222   fntExtra = new TNtuple("fntExtra","extra clusters in ITS","ncls:layer:ladder:volid:phi:x:y:z:xloc:zloc:dxy:dz:d0mu:z0mu");
 223   fListOfHistos->Add(fntExtra);
 224
 225   fntCosmicMatching = new TNtuple("fntCosmicMatching","cosmic tracks matching in ITS","ncls1:ncls2:pt1:pt2:sigmad01:sigmad02:sigmaz01:sigmaz02:dxy:dz:phimu:thetamu:d0mu:z0mu");
 226   fListOfHistos->Add(fntCosmicMatching);
 227
 228   fspTree = new TTree("spTree","Tree with ITS track points");
 229   const AliTrackPointArray *array = 0;
 230   Float_t curv,curverr;
 231   fspTree->Branch("SP","AliTrackPointArray",&array);
 232   fspTree->Branch("curv",&curv);
 233   fspTree->Branch("curverr",&curverr);
 234
 235   return;
 236 }
 237
 238 //________________________________________________________________________
 239 void AliAlignmentDataFilterITS::Exec(Option_t */*option*/)
 240 {
 241   // Execute analysis for current event:
 242   // write ITS AliTrackPoints for selected tracks to fspTree
 243
 244   // load the geometry
 245   if(!gGeoManager) {
 246     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 247     if(!gGeoManager) {
 248       printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n");
 249       return;
 250     }
 251   }
 252
 253   // check if we have AliITSRecoParam
 254   if(!AliITSReconstructor::GetRecoParam()) {
 255     if(!fITSRecoParam) {
 256       printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n");
 257       return;
 258     }
 259   }
 260
 261   if(!fESD) {
 262     printf("AliAlignmentDataFilterITS::Exec(): no ESD \n");
 263     return;
 264   }
 265   if(!fESDfriend) {
 266     printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n");
 267     return;
 268   }
 269   // attach ESDfriend
 270   fESD->SetESDfriend(fESDfriend);
 271
 272
 273   // Process event as Cosmic or Collision
 274   //if(esd->GetEventType()== ???? ) {
 275   printf("AliAlignmentDataFilterITS::Exec(): MOVE ASAP TO esd->GetEventType() !\n");
 276   if(GetRecoParam()->GetAlignFilterCosmics()) {
 277     FilterCosmic(fESD);
 278   } else {
 279     FilterCollision(fESD);
 280   }
 281
 282   return;
 283 }
 284
 285 //________________________________________________________________________
 286 void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd)
 287 {
 288   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 289   // of top and bottom track, merge the two tracks, if requested)
 290   //
 291
 292   // Set branch addresses for space points tree
 293   AliTrackPointArray *arrayForTree=0;
 294   Float_t curv,curverr;
 295   fspTree->SetBranchAddress("SP",&arrayForTree);
 296   fspTree->SetBranchAddress("curv",&curv);
 297   fspTree->SetBranchAddress("curverr",&curverr);
 298
 299   TString triggeredClass = fESD->GetFiredTriggerClasses();
 300   if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return;
 301
 302
 303   Int_t ntracks = esd->GetNumberOfTracks();
 304   if(ntracks<2) return;
 305
 306   if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return;
 307
 308   Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
 309
 310   Int_t *goodtracksArray = new Int_t[ntracks];
 311   Float_t *phiArray = new Float_t[ntracks];
 312   Float_t *thetaArray = new Float_t[ntracks];
 313   Int_t *nclsArray = new Int_t[ntracks];
 314   Int_t ngt=0;
 315   Int_t itrack=0;
 316   for (itrack=0; itrack < ntracks; itrack++) {
 317     AliESDtrack *track = esd->GetTrack(itrack);
 318     if (!track) continue;
 319
 320
 321     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 322
 323     if(GetRecoParam()->GetAlignFilterOnlyITSSATracks() && track->GetNcls(1)>0) continue;
 324     if(GetRecoParam()->GetAlignFilterOnlyITSTPCTracks() && track->GetNcls(1)==0) continue;
 325
 326     Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
 327     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
 328
 329     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 330        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 331
 332     goodtracksArray[ngt] = itrack;
 333     phiArray[ngt]        = phi;
 334     thetaArray[ngt]      = theta;
 335     nclsArray[ngt]       = track->GetNcls(0);
 336     ngt++;
 337   }
 338
 339   if(ngt<2) {
 340     delete [] goodtracksArray; goodtracksArray=0;
 341     delete [] phiArray; phiArray=0;
 342     delete [] thetaArray; thetaArray=0;
 343     delete [] nclsArray; nclsArray=0;
 344     return;
 345   }
 346
 347   // check matching of the two tracks from the muon
 348   Float_t min = 10000000.;
 349   Int_t maxCls = 0;
 350   Int_t good1 = -1, good2 = -1;
 351   for(Int_t itr1=0; itr1<ngt-1; itr1++) {
 352     for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
 353       Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
 354       if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 355       Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
 356       if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 357       if(nclsArray[itr1]+nclsArray[itr2] > maxCls) {
 358         maxCls = nclsArray[itr1]+nclsArray[itr2];
 359         min = deltatheta+deltaphi;
 360         good1 = goodtracksArray[itr1];
 361         good2 = goodtracksArray[itr2];
 362       } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) {
 363         if(deltatheta+deltaphi < min) {
 364           min = deltatheta+deltaphi;
 365           good1 = goodtracksArray[itr1];
 366           good2 = goodtracksArray[itr2];
 367         }
 368       }
 369     }
 370   }
 371
 372   delete [] goodtracksArray; goodtracksArray=0;
 373   delete [] phiArray; phiArray=0;
 374   delete [] thetaArray; thetaArray=0;
 375   delete [] nclsArray; nclsArray=0;
 376
 377   if(good1<0) return;
 378   AliDebug(2,"ok track matching");
 379
 380   // track1 will be the inward track (top)
 381   // track2 the outward (bottom)
 382   AliESDtrack *track1=0;
 383   AliESDtrack *track2=0;
 384   AliESDtrack *track = esd->GetTrack(good1);
 385   if(track->Py()>0) {
 386     track1 = esd->GetTrack(good1);
 387     track2 = esd->GetTrack(good2);
 388   } else {
 389     track1 = esd->GetTrack(good2);
 390     track2 = esd->GetTrack(good1);
 391   }
 392
 393   AliTrackPoint point;
 394   const AliTrackPointArray *array=0;
 395   Int_t ipt,volId,modId,layerId,lay,lad,det;
 396   Int_t jpt=0;
 397   Bool_t layerOK[6][2];
 398   Int_t nclsTrk[2]={0,0};
 399
 400   for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE;
 401
 402   for(itrack=0; itrack<2; itrack++) {
 403     if(itrack==0) {
 404       track = track1;
 405     } else {
 406       track = track2;
 407     }
 408     array = track->GetTrackPointArray();
 409     if(!array) {
 410       AliWarning("No tracks points avaialble");
 411       continue;
 412     }
 413     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 414       array->GetPoint(point,ipt);
 415       volId = point.GetVolumeID();
 416       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 417       AliDebug(2,Form("%d %d\n",ipt,layerId-1));
 418       if(point.IsExtra() &&
 419          GetRecoParam()->GetAlignFilterSkipExtra()) continue;
 420       if(layerId>6) continue;
 421       if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue;
 422       // check minAngleWrtITSModulePlanes
 423       Double_t p[3]; track->GetDirection(p);
 424       TVector3 pvec(p[0],p[1],p[2]);
 425       Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
 426       TVector3 normvec(rot[1],rot[4],rot[7]);
 427       Double_t angle = pvec.Angle(normvec);
 428       if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
 429       angle = 0.5*TMath::Pi()-angle;
 430       if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue;
 431       layerOK[layerId-1][itrack]=kTRUE;
 432       jpt++;
 433       nclsTrk[itrack]++;
 434     }
 435   }
 436   AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1]));
 437
 438   // read ITS cluster maps
 439   Int_t map1[6],map2[6];
 440   for(Int_t ilay=0;ilay<6;ilay++) {
 441     map1[ilay]=0; map2[ilay]=0;
 442     if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1;
 443     if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1;
 444   }
 445   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()));
 446   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()));
 447   Int_t idx1[12],idx2[12];
 448   track1->GetITSclusters(idx1);
 449   track2->GetITSclusters(idx2);
 450   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]));
 451   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]));
 452
 453
 454   if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return;
 455   AliDebug(2,Form(" Total points %d, accepted\n",jpt));
 456   fHistNpoints->Fill(jpt);
 457   fHistPt->Fill(0.5*(track1->Pt()+track2->Pt()));
 458
 459   Float_t d0z0mu[2];
 460   track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu);
 461   //printf("d0mu %f  z0mu %f\n",d0z0mu[0],d0z0mu[1]);
 462
 463   Float_t dzOverlap[2];
 464   Float_t curvArray[2],curverrArray[2];
 465   Double_t globExtra[3],locExtra[3];
 466   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks())
 467     arrayForTree = new AliTrackPointArray(jpt);
 468
 469   jpt=0;
 470   for(itrack=0; itrack<2; itrack++) {
 471     if(itrack==0) {
 472       track = track1;
 473     } else {
 474       track = track2;
 475     }
 476     curvArray[itrack] = track->GetC(esd->GetMagneticField());
 477     curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 478
 479     if(!GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
 480       jpt=0;
 481       arrayForTree = new AliTrackPointArray(nclsTrk[itrack]);
 482     }
 483     array = track->GetTrackPointArray();
 484     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 485       array->GetPoint(point,ipt);
 486       volId = point.GetVolumeID();
 487       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 488       if(layerId>6 || !layerOK[layerId-1][itrack]) continue;
 489       arrayForTree->AddPoint(jpt,&point);
 490       jpt++;
 491       switch(layerId) {
 492       case 1:
 493         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 494         break;
 495       case 2:
 496         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 497         break;
 498       case 3:
 499         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 500         break;
 501       case 4:
 502         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 503         break;
 504       case 5:
 505         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 506         break;
 507       case 6:
 508         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 509         break;
 510       }
 511       // Post the data for slot 0
 512       if(jpt==1) PostData(1,fListOfHistos); // only if this is the first points
 513       if(!point.IsExtra() ||
 514          !GetRecoParam()->GetAlignFilterFillQANtuples()) continue;
 515       nclsTrk[itrack]--;
 516       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 517       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 518       globExtra[0]=point.GetX();
 519       globExtra[1]=point.GetY();
 520       globExtra[2]=point.GetZ();
 521       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 522       //printf("%d %d %d %d %d  %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]);
 523       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 524       AliTrackPoint pointT;
 525       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 526       for(Int_t lll=0;lll<ipt;lll++) {
 527         array->GetPoint(pointT,lll);
 528         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 529         if(layerIdT!=layerId) continue;
 530         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 531         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 532         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 533         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 534         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 535         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 536         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 537         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 538         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 539         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 540         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]);
 541       }
 542     }
 543
 544     if(!GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
 545       curv = curvArray[itrack];
 546       curverr = curverrArray[itrack];
 547       fspTree->Fill();
 548     }
 549   }
 550
 551   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
 552     curv = 0.5*(curvArray[0]+curvArray[1]);
 553     curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]);
 554     /*AliTrackPoint pppt;
 555     for(Int_t iii=0;iii<arrayForTree->GetNPoints();iii++) {
 556       arrayForTree->GetPoint(pppt,iii);
 557       pppt.Dump();
 558       }*/
 559     fspTree->Fill();
 560   }
 561   PostData(0,fspTree);
 562
 563   if(!GetRecoParam()->GetAlignFilterFillQANtuples()) return;
 564   // fill ntuple with track-to-track matching
 565   Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty;
 566   Float_t d0[2],z0[2];
 567   Float_t sigmad0[2],sigmaz0[2];
 568   phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp());
 569   thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl());
 570   phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp());
 571   thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl());
 572   rotymu = TMath::ATan2(track1->Px(),track1->Pz());
 573   rotyout = TMath::ATan2(track2->Px(),track2->Pz());
 574
 575   dphi = phimu - (phiout+TMath::Pi());
 576   dtheta = thetamu - (TMath::Pi()-thetaout);
 577   if(rotymu>0) {
 578     droty = rotymu - (rotyout+TMath::Pi());
 579   } else {
 580     droty = rotymu - (rotyout-TMath::Pi());
 581   }
 582
 583   Double_t alpha = TMath::ATan2(track1->Py(),track1->Px());
 584
 585   track1->Propagate(alpha,0.,esd->GetMagneticField());
 586   track2->Propagate(alpha,0.,esd->GetMagneticField());
 587   d0[0] = track1->GetY();
 588   z0[0] = track1->GetZ();
 589   d0[1] = track2->GetY();
 590   z0[1] = track2->GetZ();
 591   Float_t dxy = -(d0[0]-d0[1]);
 592   Float_t dz  = z0[0]-z0[1];
 593   sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2());
 594   sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2());
 595   sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2());
 596   sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2());
 597   /*
 598   Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3];
 599   track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0);
 600   track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1);
 601   track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0);
 602   track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1);
 603   Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]);
 604   Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]);
 605   Float_t dxaty0 = x1aty0-x2aty0;
 606   */
 607   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]);
 608
 609   return;
 610 }
 611
 612 //________________________________________________________________________
 613 void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd)
 614 {
 615   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 616   // of top and bottom track, merge the two tracks, if requested)
 617   //
 618
 619   // Set branch addresses for space points tree
 620   AliTrackPointArray *arrayForTree=0;
 621   Float_t curv,curverr;
 622   fspTree->SetBranchAddress("SP",&arrayForTree);
 623   fspTree->SetBranchAddress("curv",&curv);
 624   fspTree->SetBranchAddress("curverr",&curverr);
 625
 626   Int_t ntracks = esd->GetNumberOfTracks();
 627
 628   if(ntracks==0) return;
 629
 630   if(esd->GetPrimaryVertexTracks()->GetNContributors()<=0) return;
 631
 632   Double_t vtxpos[3]; esd->GetPrimaryVertexTracks()->GetXYZ(vtxpos);
 633
 634   Int_t ncls=0;
 635   Double_t pt=-10000.;
 636   Double_t d0z0[2],covd0z0[3];
 637   const AliTrackPointArray *array = 0;
 638
 639   for (Int_t itrack=0; itrack < ntracks; itrack++) {
 640     AliESDtrack * track = esd->GetTrack(itrack);
 641     if (!track) continue;
 642
 643     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 644
 645     if(GetRecoParam()->GetAlignFilterOnlyITSSATracks() && track->GetNcls(1)>0) continue;
 646     if(GetRecoParam()->GetAlignFilterOnlyITSTPCTracks() && track->GetNcls(1)==0) continue;
 647
 648     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 649        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 650
 651     pt = track->Pt();
 652     ncls = track->GetNcls(0);
 653     Double_t maxd=10000.;
 654     track->PropagateToDCA(esd->GetPrimaryVertex(),esd->GetMagneticField(),maxd,d0z0,covd0z0);
 655
 656     // read ITS cluster map
 657     Int_t map[6];
 658     for(Int_t ilay=0;ilay<6;ilay++) {
 659       map[ilay]=0;
 660       if(track->HasPointOnITSLayer(ilay)) map[ilay]=1;
 661     }
 662     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()));
 663     Int_t idx[12];
 664     track->GetITSclusters(idx);
 665     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]));
 666
 667
 668     AliTrackPoint point;
 669     Int_t ipt,volId,modId,layerId,lay,lad,det;
 670     Int_t jpt=0;
 671     Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE;
 672
 673     array = track->GetTrackPointArray();
 674     if(!array) continue;
 675     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 676       array->GetPoint(point,ipt);
 677       volId = point.GetVolumeID();
 678       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 679       if(layerId<1 || layerId>6) continue;
 680       if(point.IsExtra() &&
 681          GetRecoParam()->GetAlignFilterSkipExtra()) continue;
 682       layerOK[layerId-1]=kTRUE;
 683       jpt++;
 684     }
 685
 686     if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 687
 688     fHistNpoints->Fill(jpt);
 689     fHistPt->Fill(pt);
 690     PostData(1,fListOfHistos);
 691
 692     Float_t dzOverlap[2];
 693     Double_t globExtra[3],locExtra[3];
 694     arrayForTree = new AliTrackPointArray(jpt);
 695     jpt=0;
 696     array = track->GetTrackPointArray();
 697     if(!array) continue;
 698     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 699       array->GetPoint(point,ipt);
 700       volId = point.GetVolumeID();
 701       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 702       if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue;
 703       if(!point.IsExtra() ||
 704          !GetRecoParam()->GetAlignFilterFillQANtuples()) continue;
 705       ncls--;
 706       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 707       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 708       globExtra[0]=point.GetX();
 709       globExtra[1]=point.GetY();
 710       globExtra[2]=point.GetZ();
 711       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 712       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 713       AliTrackPoint pointT;
 714       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 715       for(Int_t lll=0;lll<ipt;lll++) {
 716         array->GetPoint(pointT,lll);
 717         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 718         if(layerIdT!=layerId) continue;
 719         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 720         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 721         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 722         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 723         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 724         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 725         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 726         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 727         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 728         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 729         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]);
 730       }
 731       arrayForTree->AddPoint(jpt,&point);
 732       switch(layerId) {
 733       case 1:
 734         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 735         break;
 736       case 2:
 737         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 738         break;
 739       case 3:
 740         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 741         break;
 742       case 4:
 743         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 744         break;
 745       case 5:
 746         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 747         break;
 748       case 6:
 749         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 750         break;
 751       }
 752       jpt++;
 753     }
 754
 755     curv = track->GetC(esd->GetMagneticField());
 756     curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 757
 758     fspTree->Fill();
 759
 760   } // end of tracks loop
 761
 762   PostData(0,fspTree);
 763
 764   return;
 765 }
 766
 767 //________________________________________________________________________
 768 void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/)
 769 {
 770   // Terminate analysis
 771   //
 772   AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n");
 773
 774   fspTree = dynamic_cast<TTree*> (GetOutputData(0));
 775   if (!fspTree) {
 776     printf("ERROR: fspTree not available\n");
 777     return;
 778   }
 779
 780   fListOfHistos = dynamic_cast<TList*> (GetOutputData(1));
 781   if (!fListOfHistos) {
 782     printf("ERROR: fListOfHistos not available\n");
 783     return;
 784   }
 785
 786   fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints"));
 787   fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt"));
 788   fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0"));
 789   fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1"));
 790   fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2"));
 791   fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3"));
 792   fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4"));
 793   fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5"));
 794   fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra"));
 795   fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching"));
 796
 797
 798
 799   return;
 800 }
 801