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
  20   if (gAlice) {
  21     delete gAlice->GetRunLoader();
  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   AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
  70   if (!dITS) {
  71     cerr<<"Cannot find the ITS detector !"<<endl;
  72     return;
  73   }
  74   AliITSgeom *geom = dITS->GetITSgeom();
  75
  76   // Instance of the ITS tracker
  77   AliITStrackerSA itsTracker(geom);
  78   Int_t ITSclusters[6] = {1,1,1,1,1,1};
  79   itsTracker.SetLayersNotToSkip(ITSclusters);
  80
  81
  82   // Primary vertex reconstruction in pp
  83   AliESDVertex *initVertex = 0;
  84   TFile *invtx = new TFile("AliESDVertexMean.root");
  85   if(collcode==1 && useMeanVtx) {
  86     // open the mean vertex
  87     initVertex = (AliESDVertex*)invtx->Get("vtxmean");
  88   } else {
  89     Double_t pos[3]={0.5,0.5,0.};
  90     Double_t err[3]={3.,3.,5.};
  91     initVertex = new AliESDVertex(pos,err);
  92   }
  93   invtx->Close();
  94   delete invtx;
  95   AliVertexerTracks *vertexer = new AliVertexerTracks();
  96   vertexer->SetVtxStart(initVertex);
  97   vertexer->SetDebug(0);
  98   delete initVertex;
  99   initVertex=0;
 100
 101   /***** The TREE for ESD is born *****/
 102   TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
 103   AliESD *event=0; AliESD *eventTPCin=0;
 104   esdTree->Branch("ESD","AliESD",&event);
 105
 106   if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
 107   if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
 108   cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
 109
 110   TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
 111   AliESDVertex *vertexSPD = new AliESDVertex();
 112   Char_t zver[100];
 113   Double_t vtx[3]={0,0,0};
 114   Double_t sigmavtx[3]={0.07,0.07,0.1};
 115
 116
 117   //<---------------------------------- The Loop over events begins
 118   TStopwatch timer;
 119   Int_t trc;
 120   for(Int_t i=firstEvent; i<=lastEvent; i++) {
 121
 122     cout<<" Processing event number : "<<i<<endl;
 123     AliESD *event = new AliESD();
 124     event->SetRunNumber(gAlice->GetRunNumber());
 125     event->SetEventNumber(i);
 126     event->SetMagneticField(gAlice->Field()->SolenoidField());
 127     rl->GetEvent(i);
 128
 129     //***** Primary vertex from SPD from file
 130     sprintf(zver,"Event%d/Vertex",i);
 131     vertexSPD = (AliESDVertex*)ppZ->Get(zver);
 132     if(!vertexSPD) {
 133       esdTree->Fill(); delete event;
 134       continue;
 135     }
 136     event->SetVertex(vertexSPD);
 137     vertexSPD->GetXYZ(vtx);
 138     vertexSPD->GetSigmaXYZ(sigmavtx);
 139
 140     //***** TPC tracking
 141     if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
 142       printf("exiting TPC tracker with code %d in event %d\n",trc,i);
 143       esdTree->Fill(); delete event;
 144       continue;
 145     }
 146
 147     // make a copy of the ESD at this stage
 148     eventTPCin = event;
 149
 150     //***** ITS tracking
 151     itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
 152     TTree *itsTree=itsl->TreeR();
 153     if (!itsTree) {
 154       cerr<<"Can't get the ITS cluster tree !\n";
 155       esdTree->Fill(); delete event;
 156       return;
 157     }
 158     itsTracker.UnloadClusters();
 159     itsTracker.LoadClusters(itsTree);
 160     if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
 161       printf("exiting ITS tracker with code %d in event %d\n",trc,i);
 162       esdTree->Fill(); delete event;
 163       continue;
 164     }
 165
 166     // Bring kTPCin-tracks back to the TPC inner wall
 167     BackToTPCInnerWall(event,eventTPCin);
 168
 169     // refit inward in ITS:
 170     // - refit without vertex constraint
 171     // - propagate through beam pipe to local x = 0
 172     itsTracker.RefitInward(event);
 173
 174
 175     //***** Vertex from ESD tracks
 176     if(collcode==1) { // pp
 177       AliESDVertex *vertexTrks =
 178         (AliESDVertex*)vertexer->FindPrimaryVertex(event);
 179       event->SetPrimaryVertex(vertexTrks);
 180     }
 181
 182     esdTree->Fill();
 183     delete event;
 184
 185   }//<-----------------------------------The Loop over events ends here
 186   timer.Stop(); timer.Print();
 187
 188   // The AliESDs.root is born
 189   TFile *ef = TFile::Open("AliESDs.root","RECREATE");
 190   if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
 191
 192   //Write the tree and close everything
 193   esdTree->Write();
 194   delete esdTree;
 195   ef->Close();
 196
 197   delete vertexer;
 198   delete rl;
 199
 200   return;
 201 }
 202 //--------------------------------------------------------------------------
 203 void BackToTPCInnerWall(AliESD *event,AliESD *eventTPC) {
 204
 205   Int_t ntracks = eventTPC->GetNumberOfTracks();
 206   AliESDtrack *esdTrackTPC = 0;
 207
 208   // create relation between event and eventTPC
 209   Int_t labelsTPC[100000000];
 210   for(Int_t tr = 0; tr<ntracks; tr++) {
 211     esdTrackTPC = (AliESDtrack*)event->GetTrack(tr);
 212     labelsTPC[TMath::Abs(esdTrackTPC->GetLabel())] = tr;
 213   }
 214
 215   ntracks = event->GetNumberOfTracks();
 216   AliESDtrack *esdTrack = 0;
 217   esdTrackTPC = 0;
 218   Int_t indexTPC;
 219
 220   // loop on tracks
 221   for(tr = 0; tr<ntracks; tr++) {
 222     esdTrack = (AliESDtrack*)event->GetTrack(tr);
 223     // set to kITSout the tracks that don't have kTPCin
 224     // (they've been found by AliITStrackerSA)
 225     if(!(esdTrack->GetStatus()&AliESDtrack::kTPCin)) {
 226       esdTrack->SetStatus(AliESDtrack::kITSout);
 227       continue;
 228     }
 229
 230     // skip tracks that don't have kITSin
 231     if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
 232
 233     indexTPC = labelsTPC[TMath::Abs(esdTrack->GetLabel())];
 234     esdTrackTPC = (AliESDtrack*)eventTPC->GetTrack(indexTPC);
 235
 236     AliITStrackMI *itsTrack = 0;
 237     try {
 238       itsTrack = new AliITStrackMI(*esdTrackTPC);
 239       esdTrack = 0;
 240     }
 241     catch (const Char_t *msg) {
 242         Warning("ToTPCInnerWall",msg);
 243         continue;
 244     }
 245     itsTrack->UpdateESDtrack(AliESDtrack::kITSout);
 246     esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
 247
 248     delete itsTrack;
 249   } // end loop on tracks
 250
 251   return;
 252 }
 253 //--------------------------------------------------------------------------