TPC/AliBarrelRec_TPCparam.C

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