TPC/AliBarrelRec_TPCparam.C

   1 void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
   2   //
   3   // Macro to create AliESDs.root using parametrized TPC tracking
   4   // and AliITStrackerSA (MI + SA)
   5   //
   6   // Input files:
   7   // - galice.root
   8   // - Kinematics.root
   9   // - TrackReferences.root
  10   // - ITS.RecPoints.root (use AliRecontruction class)
  11   // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
  12   //
  13   // A. Dainese - INFN Legnaro
  14   //
  15
  16   Int_t  collcode = 1; // pp collisions
  17   Bool_t useMeanVtx = kFALSE;
  18
  19   TGeoManager::Import("geometry.root");
  20
  21   if (gAlice) {
  22     delete gAlice->GetRunLoader();
  23     delete gAlice;
  24     gAlice=0;
  25   }
  26   AliRunLoader *rl = AliRunLoader::Open("galice.root");
  27   if (rl == 0x0) {
  28     cerr<<"Can not open session"<<endl;
  29     return;
  30   }
  31   Int_t retval = rl->LoadgAlice();
  32   if (retval) {
  33     cerr<<"LoadgAlice returned error"<<endl;
  34     delete rl;
  35     return;
  36   }
  37   retval = rl->LoadHeader();
  38   if (retval) {
  39     cerr<<"LoadHeader returned error"<<endl;
  40     delete rl;
  41     return;
  42   }
  43   gAlice=rl->GetAliRun();
  44
  45
  46
  47   // Get field from galice.root
  48   AliMagF *fiel = (AliMagF*)gAlice->Field();
  49   Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
  50   // Set the conversion constant between curvature and Pt
  51   AliTracker::SetFieldMap(fiel,kTRUE);
  52
  53   /**** The TPC corner ********************/
  54
  55   AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
  56   tpcTrackerPar.Init();
  57
  58   //**** Switch on the PID class (mandatory!)
  59   AliPID pid;
  60
  61   /**** The ITS corner ********************/
  62
  63   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
  64   if (itsl == 0x0) {
  65     cerr<<"Cannot get the ITS loader"<<endl;
  66     return;
  67   }
  68   itsl->LoadRecPoints("read");
  69
  70   AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
  71   if (!dITS) {
  72     cerr<<"Cannot find the ITS detector !"<<endl;
  73     return;
  74   }
  75   AliITSgeom *geom = dITS->GetITSgeom();
  76
  77   // Instance of the ITS tracker
  78   AliITStrackerSA itsTracker(geom);
  79   Int_t ITSclusters[6] = {1,1,1,1,1,1};
  80   itsTracker.SetLayersNotToSkip(ITSclusters);
  81
  82
  83   // Primary vertex reconstruction in pp
  84   AliESDVertex *initVertex = 0;
  85   TFile *invtx = new TFile("AliESDVertexMean.root");
  86   if(collcode==1 && useMeanVtx) {
  87     // open the mean vertex
  88     initVertex = (AliESDVertex*)invtx->Get("vtxmean");
  89   } else {
  90     Double_t pos[3]={0.5,0.5,0.};
  91     Double_t err[3]={3.,3.,5.};
  92     initVertex = new AliESDVertex(pos,err);
  93   }
  94   invtx->Close();
  95   delete invtx;
  96   AliVertexerTracks *vertexer = new AliVertexerTracks();
  97   vertexer->SetVtxStart(initVertex);
  98   vertexer->SetDebug(0);
  99   delete initVertex;
 100   initVertex=0;
 101
 102   /***** The TREE for ESD is born *****/
 103   TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
 104   AliESD *event=0; AliESD *eventTPCin=0;
 105   esdTree->Branch("ESD","AliESD",&event);
 106
 107   if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
 108   if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
 109   cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
 110
 111   TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
 112   AliESDVertex *vertexSPD = new AliESDVertex();
 113   Char_t zver[100];
 114   Double_t vtx[3]={0,0,0};
 115   Double_t sigmavtx[3]={0.07,0.07,0.1};
 116
 117
 118   //<---------------------------------- The Loop over events begins
 119   TStopwatch timer;
 120   Int_t trc;
 121   for(Int_t i=firstEvent; i<=lastEvent; i++) {
 122
 123     cout<<" Processing event number : "<<i<<endl;
 124     AliESD *event = new AliESD();
 125     event->SetRunNumber(gAlice->GetRunNumber());
 126     event->SetEventNumber(i);
 127     event->SetMagneticField(gAlice->Field()->SolenoidField());
 128     rl->GetEvent(i);
 129
 130     //***** Primary vertex from SPD from file
 131     sprintf(zver,"Event%d/Vertex",i);
 132     vertexSPD = (AliESDVertex*)ppZ->Get(zver);
 133     if(!vertexSPD) {
 134       esdTree->Fill(); delete event;
 135       continue;
 136     }
 137     event->SetVertex(vertexSPD);
 138     vertexSPD->GetXYZ(vtx);
 139     vertexSPD->GetSigmaXYZ(sigmavtx);
 140
 141     //***** TPC tracking
 142     if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
 143       printf("exiting TPC tracker with code %d in event %d\n",trc,i);
 144       esdTree->Fill(); delete event;
 145       continue;
 146     }
 147
 148     // make a copy of the ESD at this stage
 149     eventTPCin = event;
 150
 151     //***** ITS tracking
 152     itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
 153     TTree *itsTree=itsl->TreeR();
 154     if (!itsTree) {
 155       cerr<<"Can't get the ITS cluster tree !\n";
 156       esdTree->Fill(); delete event;
 157       return;
 158     }
 159     itsTracker.UnloadClusters();
 160     itsTracker.LoadClusters(itsTree);
 161     if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
 162       printf("exiting ITS tracker with code %d in event %d\n",trc,i);
 163       esdTree->Fill(); delete event;
 164       continue;
 165     }
 166
 167     // Bring kTPCin-tracks back to the TPC inner wall
 168     BackToTPCInnerWall(event,eventTPCin);
 169
 170     // refit inward in ITS:
 171     // - refit without vertex constraint
 172     // - propagate through beam pipe to local x = 0
 173     itsTracker.RefitInward(event);
 174
 175
 176     //***** Vertex from ESD tracks
 177     if(collcode==1) { // pp
 178       AliESDVertex *vertexTrks =
 179         (AliESDVertex*)vertexer->FindPrimaryVertex(event);
 180       event->SetPrimaryVertex(vertexTrks);
 181     }
 182
 183     esdTree->Fill();
 184     delete event;
 185
 186   }//<-----------------------------------The Loop over events ends here
 187   timer.Stop(); timer.Print();
 188
 189   // The AliESDs.root is born
 190   TFile *ef = TFile::Open("AliESDs.root","RECREATE");
 191   if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
 192
 193   //Write the tree and close everything
 194   esdTree->Write();
 195   delete esdTree;
 196   ef->Close();
 197
 198   delete vertexer;
 199   delete rl;
 200
 201   return;
 202 }
 203 //--------------------------------------------------------------------------
 204 void BackToTPCInnerWall(AliESD *event,AliESD *eventTPC) {
 205
 206   Int_t ntracks = eventTPC->GetNumberOfTracks();
 207   AliESDtrack *esdTrackTPC = 0;
 208
 209   // create relation between event and eventTPC
 210   Int_t labelsTPC[100000000];
 211   for(Int_t tr = 0; tr<ntracks; tr++) {
 212     esdTrackTPC = (AliESDtrack*)event->GetTrack(tr);
 213     labelsTPC[TMath::Abs(esdTrackTPC->GetLabel())] = tr;
 214   }
 215
 216   ntracks = event->GetNumberOfTracks();
 217   AliESDtrack *esdTrack = 0;
 218   esdTrackTPC = 0;
 219   Int_t indexTPC;
 220
 221   // loop on tracks
 222   for(tr = 0; tr<ntracks; tr++) {
 223     esdTrack = (AliESDtrack*)event->GetTrack(tr);
 224     // set to kITSout the tracks that don't have kTPCin
 225     // (they've been found by AliITStrackerSA)
 226     if(!(esdTrack->GetStatus()&AliESDtrack::kTPCin)) {
 227       esdTrack->SetStatus(AliESDtrack::kITSout);
 228       continue;
 229     }
 230
 231     // skip tracks that don't have kITSin
 232     if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
 233
 234     indexTPC = labelsTPC[TMath::Abs(esdTrack->GetLabel())];
 235     esdTrackTPC = (AliESDtrack*)eventTPC->GetTrack(indexTPC);
 236
 237     AliITStrackMI *itsTrack = 0;
 238     try {
 239       itsTrack = new AliITStrackMI(*esdTrackTPC);
 240       esdTrack = 0;
 241     }
 242     catch (const Char_t *msg) {
 243         Warning("ToTPCInnerWall",msg);
 244         continue;
 245     }
 246     itsTrack->UpdateESDtrack(AliESDtrack::kITSout);
 247     esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
 248
 249     delete itsTrack;
 250   } // end loop on tracks
 251
 252   return;
 253 }
 254 //--------------------------------------------------------------------------